NAVAL  POSTGRADUATE  SCHOOL 
Monterey ,  California 


AD-A242  715 


DTK 

FVrCTF 

j\iOV2  ^ 


ARRIVAL  TIME  TRACKING  OF  PARTIALLY 
RESOLVED  ACOUSTIC  RAYS  WITH  APPLICATION  TO 
OCEAN  ACOUSTIC  TOMOGRAPHY 


Edwin  K.  Chaulk 


March,  1991 


Thesis  Advisor: 


James  H.  Miller 


Approved  for  public  release:  distribution  is 
unlimited . 


Ill  ilps 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No  0704-0188 


1a  REPORT  SECURITY  CLASSIFICATION 

Unclassified 


2a  SECURITY  CLASSIFICATION  AUTHORITY 


2b  DECLASSIFICATION /DOWNGRADING  SCHEDULE 


4  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 


1b  RESTRICTIVE  MARKINGS 


3  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlimited 


5  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


6a  NAME  OF  PERFORMING  ORGANIZATION 

Naval  Postgraduate  School 


6c  ADDRESS  (C/ty,  State,  and  ZIP  Code) 


6b  OFFICE  SYMBOL 
(If  applicable) 


7a  NAME  OF  MONITORING  ORGANIZATION 

Naval  Postgraduate  School 


7b  ADDRESS  (C/ty,  State,  and  ZIP  Code) 


Monterey.  California  93943-5000 


8a  NAME  OF  FUNDING /SPONSORING  8b  OFFICE  SYMBOL 

ORGANIZATION  (If  applicable) 


Monterey.  California  93943-5000 


9  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


10  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
element  NO 


project 

TASK 

no 

NO 

WORK  UNIT 
ACCESSION  NO 


8c  ADDRESS  fC/ty,  State,  and  ZIP  Code) 


11  TITLE  (/nt/ube  Security  OaM/t/cat/on;  arrival  TIME  TRACKING  OF  PARTIALLY  RESOLVED  ACOUSTIC 
RAYS  WITH  APPLICATION  TO  OCEAN  ACOUSTIC  TOMOGRAPHY 


12  personal  AUTHOR(S) 
Chaulk,  Edwin  K. 


13a  type  of  report 
Master's  Thesis 


16  jjj  jjjij  thesis  are  those  of  the  author  and  do  not  reflect  the  official 

policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 


13b  TIME  COVERED 
FROM  AuQ  BBTO  Mar  9 


T4  DATE  OF  REPORT  (Year.  Month.  Day) 
March  1991 


15  PACE  COUNT 

157 


COSATI  CODES 


SUBGROUP 


18  SUBJECT  TERMS  (Continue  on  reverse  tf  necessary  and  identify  by  block  number) 
Acoustics  Tomography,  Adaptive  Spectral 
Estimation,  Monterey  Bay,  Multipath  Resolution 


19  abstract  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Accurate  estimation  of  arrival  times  along  an  ocean  acoustic  ray  path  is  an  important  component  of 
ocean  acoustic  tomography.  A  straightforward  method  of  arrival  time  estimation,  ba,sed  on  locating 
the  maximum  value  of  an  interpolated  arrival,  was  used  with  limited  success  for  ainalysis  of  data 
from  the  December  1988  Monterey  Bay  Tomography  Experiment.  Close  examination  of  the  data 
revealed  multiple  closely  spaced  arrivals  of  similar  amplitude,  only  partially  resolved  in  many  returns. 
A  modification  to  the  original  tracking  algorithm  succeeded  in  improving  the  estimates  and  lead  to 
the  development  of  a  tracker  based  on  a  least  mean  squares  (LMS)  linear  predictive  filter.  A  second 
algorithm,  based  on  a  modified  recursive  least  squares  (MRLS)  solution,  allows  the  estimation  of 
dynamic  spectral  processes  at  surface  and  internal  wave  frequencies  in  the  tomography  arrivals. 


20  DISTRIBUTION /availability  OF  abstract  21  ABSTRACT  SECURITY  CLASSIFICATION 

unclassified/unlimited  □  same  as  rpt  □  dtic  USERS  Unclassified _ 


22a  NAME  OF  RESPONSIBLE  INDIVIDUAL  j22b  TELEPHONE  f/nc/ude  Ared  Code)  22c  OFFICE  SYMBOL 

Ralph  Hippenstiel  _ 14  08  646-26  33  _ _ _ _ 


DD  Form  1473,  JUN  86  Previous  editions  are  obsolete  _ SECURITY  classification  QF  this  page 

S/N  0102-LF-014-6603  Unclassified 


1 


Approved  for  public  release;  distribution  is  unlimited 


Arrival  Time  Tracking  of  Partially 
Resolved  Acoustic  Rays 
with 

Application  to 

Ocean  Acoustic  Tomography 

by 

Edwin  K.  Chaulk 

Lieutenant(N),  Canadian  Armed  Forces 
B.S.,  Memorial  University  of  Newfoundland,  1979 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  ELECTRICAL  ENGINEERING 


from  the 

NAVAL  POSTGRADUATE  SCHOOL 
March,  1991 


Department  of  Electrical  and  Computer  Engineering 


11 


ABSTRACT 


Accurate  estimation  of  arrival  times  along  an  ocean  acoustic  ray  path  is  an 
important  component  of  ocean  acoustic  tomography.  A  straightforward  method  of 
arrival  time  estimation,  based  on  locating  the  maximum  value  of  an  interpolated 
arrival,  was  used  with  limited  success  for  analysis  of  data  from  the  December  1988 
Monterey  Bay  Tomography  Experiment.  Close  examination  of  the  data  revealed 
multiple  closely  spaced  arrivals  of  similar  amplitude,  only  partially  resolved  in  many 
returns.  A  modification  to  the  original  tracking  algorithm  succeeded  in  improving 
the  estimates  and  lead  to  the  development  of  a  tracker  based  on  a  least  mean  squares 
(LMS)  linear  predictive  filter.  A  second  algorithm,  based  on  a  modified  recursive 
least  squares  (MRLS)  solution,  allows  the  estimation  of  dynamic  spectral  processes 
at  surface  and  internal  wave  frequencies  in  the  tomography  arrivals. 
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I.  INTRODUCTION 


This  thesis  work  is  a  continuation  of  research  resulting  from  the  1988  Mon¬ 
terey  Bay  Tomography  Experiment  (MBTE).  The  experiment  and  preliminary  data 
analysis  are  described  in  detail  in  [Ref.  1].  An  in  depth  treatment  of  the  basic  sig¬ 
nal  processing  algorithms  is  provided  in  [Ref.  2].  Early  acoustic  modeling  and  an 
environmental  assessment  is  presented  in  [Ref.  3]  wdth  more  in  depth  3-D  acoustic 
modeling  discussed  in  [Ref.  4].  Goals  of  the  1988  Monterey  Bay  tomography  experi¬ 
ment  included  quantifying  the  effects  of  surface  and  internal  waves  on  acoustic  signals 
designed  with  a  short  duration  maximal-  length  sequence,  transmitted  continuously 
from  an  ocean  acoustic  source.  This  introduction  summarizes  the  pertinent  results 
and  recommendations  of  the  original  research  [Ref.  1]  that  are  applicable  to  the  work 
to  follow. 

A.  THESIS  SUMMARY 

In  the  60  km  MBTE,  travel  time  fluctuations  were  found  to  be  five  to  ten  times 
greater  than  seen  in  previous  300  km  experiments  [Ref.  1].  Although  ray  paths  in 
the  MBTE  underwent  multiple  surface  interactions  while  the  longer  experiments  did 
not,  the  fluctuations  exceeded  the  predicted  levels  for  the  number  of  expected  surface 
interactions.  The  magnitude  of  the  arrival  time  spectra  at  surface  wave  frequencies 
did  not  agree  with  predictions,  but  the  frequency  and  spectral  shape  matched  closely 
those  computed  from  wave  buoy  data  also  collected  during  the  experiment.  The 
preliminary  analysis  estimated  the  surface  W’ave  spectral  characteristics  to  a  degree, 
but  useful  results  at  internal  wave  frequencies  could  not  be  obtained.  More  work  was 
necessary  to  characterize  the  frequency  and  amplitude  dynamics  of  the  surface  and 
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internal  wave  processes.  Both  effects  need  to  be  fully  understood  before  an  inverse 
mesoscale  mapping  of  th(.  circulation  in  the  Monterey  Bay  canyon  can  be  attempted. 

The  object  of  this  work  is  to  develop  signal  processing  algorithms  which  will 
enable  reasonably  accurate  estimation  of  the  spectral  content  of  the  data  in  both 
the  surface  and  internal  wave  frequency  domains.  It  is  desirable  to  produce  dynamic 
spectral  plots  (i.e.  variation  in  frequency  and  magnitude  over  time)  of  the  ocean 
processes  at  work  in  the  Monterey  Bay  canyon.  There  is  considerable  signal  processing 
involved  in  the  analysis  of  data  from  the  MBTE.  The  initial  processing,  including 
maximal-length  sequence  removal  or  matched  filtering  using  a  Hadamard  transform, 
is  described  in  [Ref.  1]  and  in  [Ref.  2].  Algorithms  developed  in  the  present  thesis 
work  are  applied  after  the  matched  filtering.  These  algorithms  are  of  a  general  nature 
and  can  be  adopted  to  process  any  time  series. 

Substantial  arrival  tracking  was  completed  prior  to  this  work,  but  it  is  of  limited 
usefulness  because  of  contamination  from  the  undetected  presence  of  partially  resolved 
arrivals  (i.e.  ray  paths  that  have  insufficient  temporal  spacing).  It  will  be  shown  that 
the  interference  of  the  closely  spaced  arrivals  is  responsible  for  the  anomalous  surface 
wave  magnitudes.  The  first  step  in  spectral  analysis  of  this  data  set,  was  to  improve 
the  arrival  tracking  algorithm  so  that  interference  effects  from  the  closely  spaced 
arrivals  were  minimized.  Normally,  multipath  interference  would  render  data  such 
as  these  useless.  Reasonable  results  were  obtained  by  utilizing  the  assumption  that 
the  received  ray  paths  were  stable  with  respect  to  each  other.  The  application  of 
an  adaptive  least  mean  squares  (LMS)  predictive  filter  in  a  mode  independent  of  the 
amplitude  fluctuations  of  the  received  signal,  yielded  a  considerable  improvement  in 
the  quality  of  the  arrival  time  tracks. 

Data  collection  was  restricted  to  six  hour  segments  as  dictated  by  the  capacity  of 
the  recording  media.  This  restriction  on  the  availability  of  large  contiguous  data  seg- 
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ments  coupled  with  the  lack  of  knowledge  of  internal  wave  processes  in  the  Monterey 
Bay  canyon,  prompted  the  development  of  an  adaptive  high  resolution  spectral  esti¬ 
mation  technique  that  could  handle  nonstationary  (i.e.  shifting  poles)  data  streams. 
Internal  waves,  if  present,  were  believed  to  move  through  the  region  cis  packets  or 
solitons  rather  than  having  a  well  defined  stationary  character  like  the  surface  wave 
components. 

This  thesis  focuses  on  two  primary  areas  to  process  data  from  the  MBTE,  arrival 
tracking  with  a  display  for  track  validation  and  dynamic  spectral  estimation  on  the 
tracks  in  the  surface  wave  and  internal  wave  frequency  domains.  The  two  algorithms 
developed,  are  discussed  along  with  preliminary  results  of  their  application  to  MBTE 
data  set. 


B.  THESIS  ORGANIZATION 

This  report  has  been  organized  in  the  following  way: 

1.  Chapter  II  describes  the  important  aspects  of  the  MBTE  with  an  overview 
of  the  signal  processing.  The  multipath  arrival  structure  is  investigated  and 
shortcomings  of  the  original  peak  tracking  algorithm  are  revealed. 

2.  Chapter  III  describes  the  implementation  of  the  LMS  peak  tracking  algorithm 
aJong  with  the  greyscale  track  verification  plots. 

3.  Chapter  IV  discusses  the  development  of  the  nonstationary  spectral  estimation 
procedure  along  with  some  performance  results. 

4.  Chapter  V  presents  results,  with  limited  physical  interpretation,  of  the  applica¬ 
tion  of  the  methods  developed  to  the  MBTE  data  set.  Also,  some  recommen¬ 
dations  for  future  work  and  aids  for  data  interpretation  are  discussed. 

5.  Chapter  VI  concludes  the  discussions. 
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II.  THE  MONTEREY  BAY  TOMOGRAPHY 

EXPERIMENT 


A.  DESCRIPTION 

The  MBTE  was  unique  in  that  is  was  conducted  in  a  coastal  region  with  ex¬ 
tremely  complex  bathymetry.  Ray  paths  from  transmitter  to  receiver,  transition 
steeply  from  deep  canyon  water  to  shallow  continental  shelf  water  causing  multiple 
bottom/surface  ray  interactions  in  the  shelf  region.  Figure  2.1  shows  an  example  set 
of  eigenrays,  from  transmitter  to  receiver,  generated  by  a  3-D  ray  tracing  model  [Ref. 
4].  This  set  of  rays  has  been  generated  for  Station  J,  the  primary  analysis  station 
selected  because  of  favorable  received  signal  characteristics. 

Figure  2.2  depicts  the  overall  geometry  of  the  MBTF.  The  transmitter  was 
placed  on  an  unnamed  seamount  at  LAT  36°56.3'N  and  LONG  122°17.84'W.  Nine 
receivers  were  placed  at  various  locations  as  determined  by  the  2-D  modeling  of  [Ref. 
3].  The  locations  were  spread  along  the  continentcil  shelf,  in  approximately  100  meters 
of  water,  around  the  periphery  of  the  Monterey  Bay  canyon  as  shown  in  Fig  2.2. 

The  four  goals  of  the  MBTF  were: 

1.  Investigate  experimentally  the  relation  between  the  frequency-direction  spec¬ 
trum  of  surface  waves  and  the  spectra  of  travel  time  changes  in  tomography 
signals. 

2.  Investigate  the  effect  of  internal  waves  on  tomography  signals  in  a  coastal  envi¬ 
ronment. 

3.  Investigate  the  effect  of  complex  three  dimensional  bathymetry  on  long  range 
acoustic  propagation. 

4.  Test  a  real-time  shore-based  tomography  data  acquisition  system.  [Ref.  1] 
Items  1  and  2  are  addressed  directly  in  the  work  to  follow. 
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Figure  2.1:  Sample  3-D  eigenray  solution  with  complicated  bathymetry 
along  each  path  to  station  J,  14  Dec  1988,  after  Smith. 


5 


Figure  2.2:  Source  location  (A)  and  receiver  locations  (B  -  L-2)  for  the 
MBTE. 
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Tomography  requires  a  resolvable  signal  (temporal  or  spacial  resolution)  along 
an  identifiable  (adequately  modeled)  and  stable  eigenray  path,  with  sufficient  signal- 
to-noise  ratio  (SNR)  at  the  receiver  to  process  the  arrival  time  perturbations  over  a 
long  period  [Ref.  5].  Signal  design  and  receiver  locations  were  chosen  to  optimize 
these  requirements  according  to  the  2-D  modeling  of  [Ref.  3). 

A  high-Q  omni-directional  transmitter  was  excited  to  continuously  transmit  a 
31  bit,  1.9375  sec,  maximal-length  pulse  compression  sequence.  The  bits  were  created 
by  phase  modulating  a  224  Hz  carrier  frequency  at  a  16  Hz  bit  rate  according  to  the 
equation, 

s{t)  =  cos{2Trfct  -f  A/,0)  (2.1) 

where  fc  is  the  acoustic  carrier  frequency,  t  is  time.  Mi  is  the  maximal-length  sequence 
bit  value  [-1,1]  for  the  fth  bit  and  6  is  the  phase  angle.  Signal-to-noise  ratio  of  the  de¬ 
modulated  and  compressed  received  signal  is  maximized  by  setting  6  =  tan~^{VN), 
where  N  is  the  number  of  bits  in  the  maximal-length  sequence.  The  31  bit,  1.9375 
second,  sequence  length  yields  a  Nyquist  frequency  of  0.258  Hz  for  samphng  dynamic 
ocean  processes.  Upon  demodulation  and  sequence  removal,  the  resulting  pulse  com¬ 
pression  sets  the  resolution  capability  for  fully  resolved  arrivals  to  62.5  msec,  a  single 
bit  pulse  width.  Figure  2.3  shows  a  block  diagram  of  the  major  signal  processing 
steps  utilized  for  the  received  data.  The  last  two  blocks  of  the  diagram  are  addressed 
in  the  chapters  to  follow.  For  an  in  depth  discussion  of  other  blocks,  see  [Ref.  1] 
and  [Ref.  2]. 

After  Hadamard  matched  filter  processing  or  maximal-length  sequence  removal, 
arrivals  have  an  ideal  character  as  shown  in  Fig  2.4.  Predicted  travel  times  for  station 
J  geometry  were  35  to  40  seconds.  Relative  and  not  absolute  travel  times  were  mea¬ 
sured,  since  the  sequence  repeated  every  1.9375  seconds.  Arrival  time  perturbations 
were  the  desired  measurement,  thus  absolute  travel  time  was  not  important  in  the  ap- 
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Figure  2.3:  General  steps  in  signal  processing  the  received  data  for  the 
MBTE 


Figure  2.4:  Acoustic  arrival  after  demodulation  and  matched  filtering. 
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Figure  2.5:  Two  resolvable  unambiguous  acoustic  arrivals. 

plicat  ion.  I  hc  lack  of  an  absolute  reference  docs,  however,  introduce  the  i>ossibility  of 
ambiguous  measurements.  Unambiguous  resolvable  arrivals  are  presented  in  Fig  2.5. 
Any  detectable  arrivals  with  travel  time  differences  of  more  than  1.9375  seconds  are 
ambiguous  since  they  fall  into  the  next  time  interval,  as  demonstrated  in  Fig  2.6.  A 
final  possibility  exists  for  the  arrival  structure  as  demonstrated  in  Fig  2.7.  Signals 
may  not  only  be  ambiguous,  but  may  also  be  unresolvable  or  only  partially  resolvable 
in  many  of  the  traces.  This  is  an  important  point  to  note  for  later  discussions. 

Data,  after  lladamard  matched  filter  processing,  can  be  displayed  as  in  Fig  2.8. 
This  display  is  somewhat  misleading  as  each  individual  trace  in  the  plot  is  an  average 
of  16  sejiarate  traces.  The  averaging  masks  the  fine  structure.  It  is  useful  to  show 
the  central  arrival  location  but  shows  nothing  of  what  processing  schemes  have  to 
deal  with  from  sequence  to  sequence.  An  equivalent  number  of  traces  as  displayed  in 
Fig  2.8,  are  displayed  in  Fig  2.9  without  av^eraging.  It  is  very  difficult  in  this  in.stance 
to  observe  the  dominant  trends.  The  peak  structure  is  considerably  more  comph  .\ 
than  indicated  in  Fig  2.8.  An  alternative  to  these  data  displays  is  a  greyscale  display 
borrowed  from  passive  sonar  processing  and  shown  in  Fig  2.10.  Individual  trace's,  as 
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Figure  2.8:  Waterfall  display  of  received  acoustic  signals  from  station  J, 
14  Dec  88.  Each  trace  is  31  seconds  of  data  coherently  averaged  to  one 
1.9375  second  maximal-length  sequence  period. 
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TABLE  2.1:  A  TABLE  OF  PERIODICITIES  OF  INTEREST  FOR  THE 
MBTE  DATA. 


Period 

Frequency 

Description 

5T2  sec 

0.2-0.0833  Hz 

Surface  gravity  waves  from 
fully  developed  seas 

6-22  sec 

0.1667-0.04545  Hz 

Sea  swell  periods 

1-3  min 

0.01667-0.00555  Hz 

Surf  beat 

8  min  -  24  hrs 

0.002083  -  1.157  x  10"^  Hz 

Internal  waves 
and  tides 

in  Fig  2.9,  are  quantized  into  nine  levels  and  presented  as  intensity  modulated  pixels. 
This  form  of  presentation  allows  a  large  amount  of  data  to  be  displayed  on  a  page 
and  leaves  integration  to  the  eye.  The  arrival  perturbation  dynamics  are  visible  and 
correlograms,  as  they  will  be  labeled,  are  used  later  in  the  report  with  overlayed  peak 
tracking  information  to  evaluate  tracker  performance  and  track  validity.  The  term 
correlogram  was  chosen  since,  in  this  case,  the  greyscale  represents  the  output  of  a 
matched  filter  or  equivalently,  a  correlator. 

Several  ocean  periodicities  of  interest  have  been  identified  in  [Ref.  1].  It  is 
possible  for  any  combination  of  these  periods  to  be  present  in  the  tomography  data. 
Thus,  they  are  listed  and  briefly  described  in  Table  2.1. 

B.  ESTIMATION  OF  ARRIVAL  TIMES 

The  methodology  used  in  [Ref.  1]  to  estimate  travel  times,  is  somewhat  lacking 
for  this  application.  Plots,  as  in  Fig  2.8,  were  used  to  select  what  appeared  to  be 
completely  resolved  arrivals.  The  criteria  for  determination  of  suitable  arrivals  for 
processing  were  simple.  “The  arrival  must  not  disappear  (an  indication  of  an  unstable 
path)  and  it  should  not  merge  or  split  with  another  arrival  (an  indication  that  the 
ray  paths  are  not  resolved)”  [Ref.  1).  While  these  statements  are  perfectly  valid,  the 
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criteria  were  applied  to  averaged  data  displays  on  which  a  maximum  of  one  hour  of 
data  could  be  plotted.  A  quick  look  a  Fig  2.9  clearly  shows  a  more  complex  peak 
structure  not  evident  in  the  Fig  2.8,  the  averaged  display. 

The  procedure  after  arrival  selection  was  to  identify  a  mean  track  value  from 
the  averaged  waterfall  data  plots.  A  specific  number  of  points  on  either  side  of  this 
mean  value  were  chosen  to  create  a  track  window.  A  sample  track  window,  for  station 
J  arrival  B,  is  included  in  Fig  2.8.  This  is  a  much  broader  track  window  than  was 
applied  in  the  original  processing.  Because  the  original  track  windows  were  both 
narrow  and  fixed  in  position,  it  was  possible  for  long  trends  to  move  the  arrivals 
outside  the  track  window  for  periods  during  the  six  hour  data  segment.  The  purpose 
of  the  track  window  was  to  define  a  search  region  for  an  algorithm  to  select  a  constant 
measurable  feature  on  each  arrival.  Since  only  relative  travel  times  were  of  interest, 
the  position  of  this  feature,  inside  a  1.9375  second  maximal-length  sequence  period, 
was  taken  a;s  the  arrival  time  estimate.  The  peak  amplitude  position  was  a  convenient 
feature  of  choice,  because  it  was  computationally  easy  to  locate.  This  approach  has 
proven  to  be  ineffective  on  the  MBTE  data  set,  because  of  the  presence  of  partially 
resolved  arrivals  (i.e  more  than  one  peak)  in  the  track  window. 

Basic  arrival  time  uncertainty  is  defined  in  terms  of  SNR  and  signal  bandwidth 


2nBVWR 

where  B  is  the  bandwidth  of  the  transmitted  signal  and  SNR  is  the  signal  to  noise 
ratio  [Ref.  5].  For  a  10  Db  SNR  and  a  bandwidth  of  16  Hz  the  uncertainty,  CTj,  is  3.1 
msec.  This  uncertainty  is  somewhat  reduced  because  the  quadrature  demodulation 
channels  were  sampled  at  64  Hz  which  constitutes  a  four  times  oversampling  of  the 
data  stream.  The  matched  filtering  treated  the  resulting  bit  stream  as  consisting  of 
four  separate  channels.  Appropriate  interleaving  after  matched  filtering,  permitted 
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Figure  2.11:  Arrival  periods  A,  B  and  C  in  received  acoustic  signals  from 
station  J,  14  Dec  88. 

the  use  of  curve  fitting  techniques  to  match  the  peak  trends.  Interpolation,  (cubic 

splines  selected)  of  the  curve  fits,  identified  peak  positions  (i.e.,  arrival  time)  to  less 

than  a  millisecond.  The  actual  uncertainty  is  more  than  the  interpolated  resolution, 

but  less  than  values  computed  by  Eq  2.2  for  a  reasonable  SNR.  An  average  SNR  of 

10  dB  is  a  ron.servative  estimate  of  the  actual  value. 

Figure  2.11  is  an  averaged  waterfall  display  of  station  J  for  a  period  showing 

the  three  distinct  arrival  periods  of  interest.  These  arrivals  are  designated  as  A,  B 

and  C  with  A  being  the  earliest  and  C  the  latest.  Figure  2.12  shows  the  travel  time 
( 
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Figure  2.12:  Arrival  time  track  for  arrival  A  of  station  J  using  original 
technique. 

fluctuations  or  arrival  track,  when  the  original  technique,  described  above,  is  applied 
to  arrival  A  of  station  J.  The  most  striking  feature  of  this  track,  is  the  number  of 
travel  time  estimates  tliat  appear  at  the  edges  of  the  track  window.  The  amount  of 
clipping  appears  excessive  since  the  averaged  waterfall  plots  such  as  Fig  2.11  would 
indicate  that  more  of  the  estimates  should  fall  into  the  track  window. 

The  complex  peak  structure,  indicated  in  Fig  2.9,  prompted  a  check  of  the 
premise  that  only  single  peaks  exist  in  a  track  window.  The  verification  procedure 
utilizes  a  slightly  different  approach  to  the  problem.  To  accommodate  Fast  Fourier 
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Transform  (FFT)  interpolation,  the  number  of  points  in  the  track  window  was  se¬ 
lected  to  be  a  power  of  two  (i.e.  8,  16,  32,  etc.).  Interpolation  using  the  FFT  is 
accomphshed  by  Fourier  transforming  a  sequence  and  zero  padding  the  center  of  the 
result  to  the  desired  power  of  two.  Performing  an  inverse  Fourier  transform  on  the 
padded  sequence  yields  the  new  interpolated  sequence,  which  in  this  case,  yields  a 
smooth  curve  that  can  be  numerically  differentiated  with  acceptable  accuracy.  Inter¬ 
polated  track  windows  are  doubly  differentiated  to  enable  removal  of  local  minima 
since  only  the  local  maxima  or  peaks  are  of  interest.  Fig  2.13  shows  the  performance 
of  this  processing,  for  a  single  trace  from  station  J  arrival  B,  with  a  16  point  repre¬ 
sentative  track  window  interpolated  to  512  points.  The  solid  line  is  an  overlay  of  the 
interpolation  on  the  track  window  points,  which  are  connected  by  the  dashed  line. 
The  computed  positions  of  the  local  maxima  are  indicated  by  the  vertical  lines.  It  is 
important  to  note  that  although  one  peak  is  dominant,  more  than  one  is  present. 

The  problem  with  ptirtiaJly  resolved  arrivals  is  cleaxly  demonstrated  in  the  next 
two  figures.  Figure  2.14  and  Fig  2.15  show  two  consecutive  mciximal-length  sequence 
lengths  or  two  consecutive  data  frames  from  station  J.  The  track  window  is  high¬ 
lighted  by  overlaying  the  FFT  interpolation  as  demonstrated  in  Fig  2.13.  Vertical 
lines  mark  the  peak  positions  as  computed  using  the  second  derivative.  Figure  2.14 
shows  a  dominant  peak  near  the  center  of  the  track  window.  In  Fig  2.15,  the  domi¬ 
nant  peak  has  moved  to  the  left  of  the  track  window  and  is  very  pronounced.  Note 
however,  the  presence  of  a  very  distinct  low  level  peak  at  the  approximate  position 
of  the  dominant  peak  of  the  previous  figure.  The  peak  amplitude  tracker  described 
above  would  indicate  an  arrival  time  shift  between  dominant  peaks  of  the  adjacent 
data  frames. 

This  measured  shift  is  falst  Interference  effects  between  the  partially  resolved 
arrivals  cause  large  amplitude  fluctuations  in  the  arrival  structure,  making  each  of 
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Figure  2.13;  Station  J  arrival  B  16  point  track  window  FFT  interpolated 
to  512  points,  showing  the  positions  of  the  local  maxima  located  by  a 
numerical  second  derivative  operation. 


Figure  2.14:  Multiple  peak  structur 
ber  204. 


Amplitude 


the  closely  spaced  arrivals  amplitude  dominant  in  different  data  frames.  There  is  no 
doubt  that  the  amplitude  fluctuations  are  driven  by  path  differences  induced  by  such 
processes  as  the  ray  interactions  with  the  surface  wave  structure,  but,  the  arrival  time 
estimation,  as  implemented,  constitutes  a  non-linear  filter.  In  this  case,  the  surface 
wave  component  of  the  signal  is  enhanced  by  actually  measuring  shifts  of  amplitude 
dominance  between  closely  spaced  arrivals.  This  explains  the  excessive  travel  time 
fluctuations  noted  for  the  surface  wa\'e  frequencies  in  [Ref.  1). 

Figure  2.16  was  produced  for  arrival  B,  the  highest  amplitude  station  J  arrival, 
by  running  the  FFT  interpolation  and  second  derivative  routines  on  the  track  window 
for  each  frame  of  the  data  and  computing  a  histogram  from  all  the  peak  measure¬ 
ments.  The  figure  shows  at  least  seven  distinct  peaks  in  this  window.  The  peaks  on 
the  extreme  left  and  right  sides  of  the  window  might  be  dismissed  as  edge  effects  from 
the  interpolation.  This  would  still  leave  five  arrivals  which  contribute  to  the  partial 
resolution  problem.  The  arrival  fluctuations  of  interest  in  tomography  would  be  de¬ 
viations  about  a  single  peak  of  the  histogram.  Obviously,  it  is  desirable  to  develop  an 
algorithm  that  would  sort  through  the  peak  structure,  independent  of  the  amplitude 
of  the  peaks  in  the  track  window  for  each  data  frame,  and  select  the  peak  belonging 
to  the  same  sub-arrival  as  in  the  previous  frame.  This  is  the  basic  concept  used  in 
the  advanced  tracker  discussed  in  the  next  chapter. 


III.  LMS  ARRIVAL  TIME  TRACKING 


Since  the  difficulties  with  the  original  tracking  method  were  caused  by  its  depen- 
dance  on  the  absolute  amplitude  of  the  arrival  peaks,  a  tracker  that  ignores  absolute 
amplitude  information  should  give  better  performance.  The  FFT  interpolation  and 
second  derivative  method  for  locating  local  maxima,  as  described  in  earlier,  is  quite 
adequate  to  give  accurate  peak  arrival  time  estimates  of  all  the  peaks  in  a  track  win¬ 
dow  for  each  maximal-length  sequence  period.  What  is  needed,  is  a  way  to  pick  the 
correct  sub-arrival  from  trace  to  trace.  At  first,  a  simple  exponential  average  of  the 
estimated  arrival  time  was  computed  at  each  step.  The  peak  selected  in  the  next  track 
Weis  the  peak  with  the  closest  arrival  time  to  the  average  value  being  maintained  from 
the  previous  traces.  This  produced  a  simple  adaptive  scheme  that  showed  improved 
results.  The  variance  on  the  track  was  reduced  amd  so  wais  the  clipping.  An  adap¬ 
tive  predictor  algorithm  would  be  capable  of  following  more  complex  fluctuations  and 
trends  in  the  data  than  the  simple  exponential  average.  The  new  concept  utilizes  a 
Widrow-Hoff  least  mean  squares  adaptive  algorithm  with  an  efficient  implementation 
to  replace  the  exponential  average  tested  in  early  development. 

A.  OPERATION  OF  THE  LMS  PEAK  TRACKING  ALGORITHM 

The  Widrow-Hoff  least  mean  squares  adaptive  algorithm  is  discussed  in  [Ref. 
6].  Its  implementation  as  an  adaptive  line  enhancement  technique  is  investigated 
thoroughly  in  [Ref.  7].  In  the  present  application,  the  algorithm  operates  akin  to  a 
sort  of  phase  lock  loop.  The  algorithm  is  initialized  in  the  vicinity  of  a  signal  and  it 
is  required  to  sor:  out  the  dominant  process  and  adapt  its  coefficients  to  follow  this 
process  as  closely  as  possible.  The  algorithm  uses  as  a  measurement,  the  arrival  peak 
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position  closest  to  the  LMS  predicted  position.  The  adaptation  feature  allows  it  to 
follow  long  trends  in  the  data.  In  this  form  of  implementation,  the  LMS  algorithm 
will  give  preference  to  the  signal  with  the  least  dynamics. 

Figure  3.1  is  a  schematic  block  diagram  of  the  LMS  implementation  for  the 
arrival  tracking  problem.  The  algorithm  has  two  adjustable  parameters,  the  order  of 
the  tap-weight  coefficients  and  the  adaptation  parameter  y,  in  the  tap-weight  vector 
update  equation.  This  equation  is  given  as, 

W{k-\-l)  =  W{k)  -  2ye{k)X{k)  (3.1) 

where  W(k  -f  1)  is  the  tap- weight  vector  to  be  applied  in  the  A:  -f  1  iteration,  W_{k) 
is  the  present  tap-weight  vector,  y  is  the  adaptation  parameter  mentioned  above  and 
e{k)  is  the  error  between  the  prediction  and  the  measurement  for  the  kih  iteration. 
The  predictor  operates  as  a  conventional  finite  impulse  response  (FIR)  filter,  with 
the  application  of  a  tap-weight  vector  to  L,  where  L  is  the  filter  order,  previous 
measurements  of  the  process.  The  difference  is  that  a  prediction  error  filter  is  formed 
and  the  weights  are  adjusted  to  minimize,  in  a  least  mean  squares  sense,  the  error 
between  the  prediction  and  the  mccisurement.  This  type  of  filter  is  able  to  handle 
nonstationary  data  streams  with  slowly  shifting  poles. 

There  are  three  items  which  must  be  addressed  to  use  this  algorithm  effectively 
in  this  application.  The  first  is  how  to  set  the  filter  length,  the  second  is  how  to 
determine  a  value  for  y  and  finally  how  is  the  filter  to  be  initialized.  The  LMS 
algorithm  is  a  relatively  inexpensive  computation,  so  the  filter  order  can  selected 
based  on  data  characteristics.  Obviously,  the  longer  the  filter,  the  more  of  the  past 
me.asurements  that  will  be  incorporated  into  the  prediction.  Ninety-six  coefficients, 
utilizing  3.1  minutes  of  past  data,  provides  sufficient  memory  for  the  first  three  periods 
in  Table  2.1.  Ideally  a  comparison  study  of  the  effects  of  filter  length  should  be 
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Figure  3.1:  Schematic  block  diagram  of  the  Widrow-Hoff  LMS  filter  im 
plementation  for  the  arrival  tracking  application. 


completed.  This  must  be  done  when  the  data  and  algorithms  have  been  moved  to  a 
more  capable  machine.  Present  run  times  are  on  the  order  of  24  hours  for  a  single 
arrival  track.  This  is  not  much  longer  than  the  original  method  and  considering 
the  improvement  that  will  be  demonstrated  later,  the  results  warrant  the  additional 
processing  time,  but  make  a  full  optimization  inappropriate. 

The  adaptation  parameter  can  be  chosen  in  accordance  with  some  simple  rela¬ 
tions  derived  in  [Ref.  7].  The  full  derivation  of  the  the  Widrow-Hoff  LMS  algorithm, 
with  the  solution  via  the  method  of  steepest  descent,  is  not  included  in  this  thesis 
work.  Sufficient  literature  exists  on  most  aspects  of  the  implementation  and  the  filter 
characteristics,  for  a  variety  of  applications.  Chapter  V  of  [Ref.  6]  is  devoted  to  the 
development  of  adaptive  algorithms  based  on  LMS  techniques.  It  is  sufficient  here  to 
state  some  of  the  more  important  results. 

The  adaptation  time  constant  is  related  inversely  to  the  eigenvalues  of  the  cor¬ 
relation  matrix  of  the  process.  There  are  as  mamy  time  constants  as  there  are  filter 
weights  according  to. 
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‘p 


4/^Ap 


(3.2) 


where  Tp  is  the  pth  adaptation  time  constant,  p  is  the  selectable  adaptation  param¬ 
eter  and  Ap  is  the  pth  eigenvalue  of  the  correlation  matrix  [Ref.  7).  Some  further 
manipulation  will  show  that. 


4pAave  4fitr{R)  ^  ^ 

where  is  the  convergence  time  constant  of  the  mean  square  error,  Xave  is  the 
average  eigenvalue  cind  tr{R)  is  the  trace  of  the  correlation  matrix.  The  significance 
the  of  Eq  3.3,  ia  that  it  chows  that  p  must  be  chosen  less  than  IjXmax  for  the  filter 
to  converge  [Ref.  7].  The  closer  p  is  chosen  to  this  value  the  faster  the  convergence, 
but  also  the  more  misadjustment  noise  occurs  in  the  weight  vector  and  thus,  more 
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noise  is  generated  in  the  filter  output. 

For  the  MBTE  data  set  the  parameter  is  chosen  to  achieve  a  desired  result 
during  the  filter  initialization  procedure.  First  note  that  this  filter  is  required  to 
operate  with  non-zero  mean  track  values.  In  fact,  the  performance  appears  to  be 
better  with  the  mean  offset,  than  if  the  track  is  processed  with  the  mean  removed. 
This  point  is  not  clearly  understood,  but  is  related  to  the  generation  of  the  W_  tap- 
weights.  The  values  of  the  weights  generated,  are  of  greater  magnituo^  fo-  a  non-zero 
mean  offset,  than  those  generated  with  a  zero  mean  sequence.  The  larger  tap-weight 
values  tend  to  make  the  filter  more  responsive  in  this  application.  It  may  be  that 
the  zero  mean  filter  simply  needs  more  time  to  properly  adapt  to  the  sequence.  In 
any  case,  the  filter  is  initialized  by  beginning  with  a  filter  of  one  coefficient  and 
adding  a  coefficient,  to  increase  the  filter  size,  as  each  new  data  sample  is  processed. 
This  is  continued  until  the  desired  filter  order  is  met.  The  adaptation  parameter  is 
selected  such  that,  by  the  time  the  desired  filter  order  has  been  achieved,  the  start 
up  transient  has  reached  a  desired  mean  arrival  time.  The  algorithm  never  failed  to 
converge  using  this  criteria,  indicating  that  the  n  <  1/A„,ai  limitation  mentioned 
earlier  has  not  been  violated.  This  method  of  setting  ^  might  not  work  for  other 
filter  orders  or  mean  arrival  times  but  is  a  good  first  try  for  most  cases. 

During  the  initialization  process,  the  algorithm  is  forced  to  chose  values  closest 
to  a  fixed  mean  line  of  interest  that  has  been  predetermined  :id  is  one  of  the  inputs. 
As  soon  cis  the  number  of  tap- weights  reaches  the  desired  order,  the  algorithm  is  left 
to  run  on  its  own  predictions.  The  start  up  transient  can  be  eliminated  by  running 
the  filter  backward  on  the  data  after  the  process  has  proceeded  forward  for  some 
time.  In  fact,  for  better  accuracy,  the  algorithm  could  be  set  to  proceed  forward  and 
backward  through  the  data  until  some  specified  global  error  tolerance  between  passes 
is  achieved.  A  future  implementation  of  this  sort  could  provide  some  interesting 
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results.  Figure  3.2  shows  the  start  up  transient  overlayed  on  a  correlogram  for  arrival 
B  of  station  J.  About  100  points  are  lost  for  post-processing.  This  is  not  significant 
enough,  at  this  point  in  the  concept  development,  to  warrant  implementation  of  the 
backward  filter  to  eliminate  the  transient.  Production  code  for  processing  arrivals 
from  all  the  stations  of  the  MBTE,  should  include  this  feature. 

A  simple  test  case  for  the  algorithm  was  devised  by  creating  a  short  data  set 
consisting  of  three  parallel  noisy  signals,  spaced  approximately  as  in  the  station  J  data 
set.  A  fourth  noisy  sinusoid  with  an  amplitude  spanning  all  three  parallel  signals  was 
added.  The  task  of  the  LMS  algorithm  was  to  track  the  sinusoidal  test  signal  through 
the  contamination  caused  by  the  other  paths.  The  test  data  set  is  shown  in  Fig  3.3. 
The  results  of  the  tracking  is  shown  in  Fig  3.4.  The  results  are  quite  reasonable.  Some 
capture  by  each  of  the  parallel  paths  is  evident,  but  the  algorithm  succeeds  in  the 
tracking  the  sinusoid  with  some  phase  delay  as  would  be  expected  from  any  filtering 
operation.  Bearing  in  mind  that  the  noise  levels  selected  for  the  test  are  quite  severe 
as  seen  in  Fig  3.3,  distortion  could  be  minimized  by  a  second  low  peiss  filter  operation 
designed  to  smooth  the  effects  of  the  path  capture  experienced  by  the  tracker  for  the 
test  data.  Variation  of  the  parameters  of  the  tracker  could  also  improve  performance. 
The  parameters  used  in  the  test  processing  were,  a  filter  order  of  96  tap-weights  and 
an  adaptation  parameter  of  0.003. 

B.  COMPARISON  OF  ARRIVAL  TRACKING  METHODOLOGIES 

At  this  point,  it  is  useful  to  compare  arrival  tracking  techniques  to  show  the 
similarities  and  differences.  The  original  processing  from  matched  filtering  to  post 
processing  is  shown  in  the  schematic  block  diagram  of  Fig  3.5.  This  figure  emphasizes 
many  aspects  of  the  processing.  The  right  hand  side  blocks  of  the  figure  indicate 
some  of  the  storage  requirements,  the  hardware  used  and  the  software  necessary 
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Figure  3.5;  Schematic  block  diagram  of  original  peak  tracking  algorithm 
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for  the  method.  All  of  the  processing  for  the  original  method  was  done  on  a  PC- 
AT  compatible  Zenith  computer  with  interfaces  to  BernoulH  21  Mb  macs  storage 
removable  hard  disks.  Most  of  the  programs  were  written  in  Microsoft  FORTAN 
Version  4.0  for  the  PC  and  are  included  in  [Ref.  1].  MATLAB,  a  product  of  The 
Mathworks  Inc.,  Sherborn,  MA  and  SURFER,  a  product  of  Golden  Software  Inc., 
Golden,  CO  were  used  mainly  for  their  plotting  capabilities  in  the  original  processing. 
One  of  the  limitations  to  analyze  this  data  set,  was  the  restrictions  on  processing 
time,  memory  and  storage  imposed  by  the  use  of  the  PC.  The  availability  of  more 
powerful  machines  such  as  the  SUN  workstations  will  make  future  processing  much 
more  effective. 

Figure  3.6  contrasts  the  implementation  of  the  LMS  method  with  the  original 
method  of  Fig  3.5.  Note  that  a  VAX  11/785  virtual  memory  machine  is  used  for  some 
of  the  processing  steps  in  the  new  implementation.  The  need  to  manipulate  larger 
arrays  outstripped  the  capabilities  of  the  PC  and  thus  MATLAB  on  the  VAX  was 
quite  useful  in  this  respect.  Additionally,  the  program  for  generating  the  correlogram 
displays  was  part  of  a  larger  software  package  written  in  FORTRAN  for  the  VAX 
computer  with  the  only  supported  output  device  being  the  Imagen  laserprinter  con¬ 
nected  to  the  VAX  machine.  The  new  processing  programs  were  written  in  MATLAB 
and  are  directly  portable  to  any  machine  running  the  MATLAB  software  package. 
These  programs  are  simple  to  understand,  easy  to  write  and  utilize  double  preci¬ 
sion  arithmetic  at  all  stages  of  processing.  To  solve  a  data  access  problem  with  the 
Bernoulli  disks,  a  MATLAB  MEX  file  interface  program  was  written  in  FORTRAN 
to  directly  interface  the  Bernoulli  disks  with  the  MATLAB  software. 

Besides  the  differences  in  display  of  the  match  filtered  output  (correlogram  vs 
SURFER  produced  waterfalls),  the  original  method  used  a  fixed  spline  interpolation 
procedure  versus  the  FI  T  interpolation  procedure  adopted  for  the  LMS  method.  The 
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Figure  3.6;  Sclieniatic  block  diagram  of  LMS  peak  tracking  algorithm 


limitations  of  the  spline  method  implemented  were  discussed  in  Chapter  II  section  B. 
The  FFT  interpolation  in  the  LMS  implementation,  has  the  limitation  that  the  num¬ 
ber  of  points  in  the  track  window  must  be  a  power  of  two.  The  interpolated  resolution 
required  for  the  FFT  must  also  be  a  power  of  two.  Providing  these  limitations  are 
met,  the  absolute  amplitude  independence  of  the  differentiated  data  means  that  the 
track  window  can  be  large  as  desired  without  fear  of  interference  from  close  resolved 
arrivals,  as  is  the  case  with  the  peak  amplitude  method.  Experimentally,  it  hais  been 
determined  that  a  good  combination  for  this  data  set  is  a  16  point  track  window  with 
a  512  point  FFT  interpolation. 

Another  significant  difference,  besides  the  arrival  time  selection  method  which 
will  be  discussed  later,  is  the  form  of  the  output  from  the  two  different  approaches. 
There  is  more  information  available  from  the  LMS  implementation.  Two  LMS  filters 
are  run  in  parallel;  one  using  the  arrival  time  information  of  the  peaks  and  the  other 
using  the  amplitudes  of  the  peaks  measured  by  the  arrival  time  filter.  A  feature  of 
the  LMS  filter  is  that  it  automatically  divides  the  tracks  into  high  and  low  frequency 
regions.  The  LMS  filter  predictive  output,  tracks  slow  trends  in  the  data,  while  the 
difference  between  the  measured  values  and  the  predicted  values  or  the  so  called  error 
signal,  track  the  high  frequency  components  of  the  data.  In  terms  of  this  data  set, 
the  arrival  time  filter  places  the  higher  frequency  surface  wave  components  in  the 
error  signal,  while  providing  information  about  the  internal  wave  spectral  region  in 
the  predictive  filter  output.  Although  the  LMS  algorithm  operating  on  the  amplitude 
is  restricted  to  utilize  the  values  measured  by  the  arrival  time  filter,  it  can  still  be 
used  to  form  a  simple  track  quality  figure.  Separate  from  the  actual  amplitude  LMS 
filter  output  tied  to  the  arrival  time  LMS  selections,  the  routine  is  allowed  to  select 
a  peak  from  the  differentiated  peak  set.  This  selection  is  compared  to  the  arrival 
time  selection.  If  the  same  peak  is  selected,  a  counter  is  incremented.  A  count  of 
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the  number  of  times  the  filters  chose  the  same  peak  from  the  differentiated  peak 
set,  divided  by  the  toted  number  of  frames  processed,  defines  a  percentage  value  of 
rehability.  If  the  processes  are  in  fact  predictable  and  the  filter  is  in  fact  tracking,  these 
two  filters  should  select  the  same  point  a  majority  of  the  time.  A  running  performance 
indication  is  then  available  as  the  tracker  proceeds  through  the  data.  Additionally, 
phase  measurements  are  recorded  for  the  arrival  times  but  an  LMS  predictive  filter  was 
not  applied.  An  LMS  filter  operating  on  the  phase  information  is  easily  implemented, 
but  the  significance  of  the  phase  as  it  applies  to  the  physical  process  is  not  yet  well 
understood.  It  is  expected  that  phase  values,  once  well  understood,  will  play  an 
important  role  in  improving  the  tracking  ability  of  the  LMS  technique. 

The  original  method  provides  four  outputs;  arrival  time,  arrival  phase,  arrival 
amplitude  and  a  track  quality  figure  based  on  a  SNR  measurement.  This  SNR  mea¬ 
surement  is,  however,  of  little  value  considering  the  interference  of  the  multiple  peaks 
in  the  arrival  windows.  Post-processing  in  both  cases  refers  to  spectred  estimation  in 
the  surface  and  internal  frequency  domains.  As  mentioned,  data  from  the  error  signal 
of  the  LMS  technique  lends  itself  to  immediate  spectral  processing  in  the  surface  wave 
region.  The  maximum  amplitude  in  the  track  window  criterion  requires  low  pass  fil¬ 
tering  to  make  the  track  usable  in  any  region,  since  the  high  frequency  effects  of  the 
clipping  have  to  be  eliminated.  Figure  3.7  is  a  flow  diagrcim  of  the  tracking  process 
indicating  three  arrival  selection  criteria.  Figure  2.12,  Fig  3.8  and  Fig  3.9  show  the 
processing  results  for  each  criterion  as  it  is  applied  to  arrivals  of  station  J.  To  allow  a 
fair  comparison.  Fig  3.9  shows  the  peak  positions  derived  from  the  second  derivative 
operation  using  the  LMS  filter  and  not  the  output  of  the  filter  itself  which  is  shown 
in  Figure  3.10.  The  second  derivative  maximum  peak  amplitude  method  shown  in 
Fig  3.7,  shows  some  slight  improvement  when  compared  to  the  track  character  of  the 
mzocimum  amplitude  in  the  window  method  shown  in  Fig  2.12,  but  the  improvement 
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Figure  3.7:  Schematic  block  diagram  for  comparison  of  peak  selection 
methods. 
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Figure  3.8:  Arrival  B  station  J  track  generated  by  the  second  derivative 
peak  amplitude  method. 

in  Fig  3.9  is  much  more  dramatic.  Clipping  has  been  eliminated  and  the  variance  on 
the  track  is  considerably  reduced.  The  results  of  the  filter  output  of  Fig  3.10  are  even 
more  impressive. 

An  effort  was  made  to  have  the  LMS  algorithm  attempt  to  respond  to  the  dom¬ 
inant  arrival  in  the  track  window  by  biasing  the  results  toward  the  higher  amplitude 
arrivals  in  the  window.  This  biasing  was  achieved  by  forcing  the  second  derivative 
algorithm  to  ignore  peaks  below  a  level  determined  by  the  average  amplitude  values 
of  all  the  peaks  in  the  window  for  a  particular  trace.  The  amplitude  biasing  has  a 
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Figure  3.9:  Arrival  B  station  J  track  generated  by,  LMS  adaptive  filter 
selecting  the  closest  peak  in  the  differentiated  peak  set. 
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Figure  3.11:  Arrival  B  station  J  track  generated  at  the  output  of  the  LMS 
adaptive  filter  with  applied  amplitude  bias  in  the  arrival  selection. 

dramatic  effect  and  shows  the  peak  switching  phenomena  of  the  original  technique 
much  clearer.  Compare  the  arrival  B  track  in  Fig  3.11,  generated  by  the  LMS  algo¬ 
rithm  with  amplitude  bias,  with  that  of  Fig  3.10,  generated  by  the  LMS  algorithm 
without  amplitude  bias.  Figure  3.11  shows  a  periodic  switching  behavior  which  is  a 
result  of  an  interaction  between  closely  spaced  arrivals. 

The  final  check  on  performance  comes  from  processing  arrivals  A,  B  and  C  of 
station  J  and  overlaying  the  results  on  the  appropriate  correlograms.  Three  LMS 
tracks  are  overlayed  on  the  correlograms.  Figure  3.12  show's  the  initial  20  minute 
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data  block  of  station  J.  Appendix  D  contains  the  correlograms  for  all  the  data  pro¬ 
cessed  from  station  J.  Next,  tv'o  additional  tracks  are  overlayed  for  the  performance 
comparison.  These  additional  overlays  consist  of  a  low  pass  filtered  arrival  A  track 
from  the  maximum  amplitude  in  the  window  method  (see  Fig  2.12)  and  a  low  pass  fil¬ 
tered  arrival  B  track  from  the  second  derivative  peak  amplitude  method  (see  Fig  3.8). 
Although  these  traces  are  all  plotted  with  solid  lines,  they  are  easily  distinguished  by 
character  alone.  The  LMS  tracks  all  have  a  fine  structure.  The  other  two  tracks  are 
smooth  with  considerable  more  variance  than  the  LMS  tracks.  The  improvements  the 
LMS  technique  provides  are  evident  in  the  correlograms  of  Appendix  D.  The  LMS 
technique  tracks  the  fine  structure  much  more  closely  than  the  other  methods.  It  is 
also  evident  that  the  method  performs  much  better  on  the  highest  amplitude  arrivals. 
The  arrival  B  track  is  much  better  behaved  than  the  arrival  A  of  C  tracks  which  have 
periods  of  low  SNR  or  total  lack  of  signal.  The  relationship  of  the  track  to  the  data 
is  quite  clear  in  the  correlogram  plots  and  as  such,  they  provide  a  necessary  check 
on  track  validity  which  will  be  required  in  the  interpretation  phase  of  data  from  the 
MBTE. 

C.  LMS  TRACKING  SUMMARY 

The  advantages  of  the  LMS  tracking  algorithm  can  be  summarized  as  follows: 

1.  No  dependence  on  the  absolute  amplitude  peak  of  the  arrival  structure. 

2.  Learns  about  the  process  eis  it  proceeds,  thus  the  track  quality  improves  with 
each  update. 

3.  The  implementation  is  extremely  fast  and  simple. 

4.  Automatically  divides  the  data  into  frequency  domains,  in  this  case,  surface 
and  internal  wave  regions. 

5.  Provides  adjustable  parameters  to  (filter  length  and  adaptation  parameter)  to 
adjust  tracking  as  required  by  the  arrival  structure. 

6.  Does  not  require  a  fixed  track  window  after  initialization,  and  thus,  can  adapt 
to  signals  that  vary  slowly  over  a  wide  region. 
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The  one  disadvantage  for  this  tracker  is  that  there  is  no  easy  way  to  control 
which  sub-arrival  the  routine  locks  on  to.  One  must  be  satisfied  with  the  routine  s 
choice  or  set  up  an  iterative  forward-backward  implementation  with  multiple  passes 
to  achieve  a  desired  result. 
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IV.  MODIFIED  RLS  SPECTRAL  ESTIMATION 


Wave  buoy  data,  collected  during  the  Monterey  Bay  experiment,  provide  a  mea¬ 
surement,  against  which,  spectral  processing  of  the  tomography  data  can  be  compared 
for  the  surface  wave  frequency  domain.  No  such  truth  exists  for  any  internal  wave 
phenomena  that  might  have  also  been  present  during  the  experiment.  To  handle  the 
possibility  that  internal  wave  mecisurements  might  be  highly  nonstationary,  an  adap¬ 
tive  algorithm  was  devised  to  provide  snapshots  of  the  process  at  a  resolution  better 
than  normally  available  from  conventional  periodogram  based  processing.  The  tech¬ 
nique  is  a  combination  of  two  algorithms  described  in  [Ref.  6].  A  modified  forward- 
backward  linear  predictor  (MFBLP)  is  implemented  using  an  update  methodology 
borrowed  from  an  recursive  least  squares  (RLS)  technique.  An  adjustable  forgetting 
factor  enables  the  algorithm  to  handle  both  stationary  or  nonstationary  (shifting 
poles)  data  streams.  This  algorithm  is  discussed  along  with  some  performance  test 
results  in  the  following  sections.  The  combined  algorithm  provides  identical  results 
as  the  MFBLP  algorithm  [Ref.  6]  for  stationary  fixed  length  data  sequences.  For  an 
in  depth  comparison  of  the  MFBLP  technique  with  other  modern  spectral  estimation 
techniques  see  [Ref.  8]. 

A.  METHOD  OF  LEAST  SQUARES 

In  order  to  combine  algorithms,  some  commonality  must  exist.  The  MFBLP  and 
RLS  algorithms  are  connected  by  a  common  basic  equation  requiring  a  least  squsires 
solution.  The  difference  is  in  how  the  least  squares  solution  is  obtained  in  each  case. 
The  RLS  method  depends  on  the  matrix  inversion  lemma  [Ref.  6,  page  385]  which 
yields  a  recursive  implementation,  while  the  MFBLP  depends  on  a  minimum  norm 
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solution  ^^roduced  from  a  matrix  pseudoinverse  [Ref.  6,  pages  336-337].  In  earh  case, 
a  tap- weight  vector  for  an  autoregressive  process  is  produced.  This  tap  weigl.l  vector 
is  applied  to  the  sample  space  to  form  a  linear  predictor.  This  predictor  is  used  to 
extend  a  data  subset  in  a  section  of  interest.  Forming  a  periodogram  of  the  extended 
data  yields  a  high  resolution  snapshot  of  the  process  for  the  selected  data  section. 
Collecting  and  displaying  these  snapshots,  at  some  interval  in  a  waterfall  or  sonogram 
display  (sonogram  refers  to  frequency  greyscale  displays  while  correlogram  refers  to 
the  output  of  a  matched  filter  in  greyscale  format),  produces  a  picture  of  the  spectral 
dynamics  of  the  process. 

The  deterministic  normal  equation  for  the  least  squares  problem  is  given  by, 


Aw  =  A^b.  (4.1) 

where  H  is  the  Hermitian  operator  or  complex  conjugate  transpose,  A  is  the  data 

mo.,,*x,  the  tap- weight  vector  for  which  the  sum  of  error  squares  is  a  . n, 

A^  is  a  forward- backward  data  matrix  as  given  by, 

\  x{M)  •••  x{N-\)  \  x*(2)  •••  x-(iV-M  +  l)l 

^  I  x(A/-l)  ...  x{N-2)  \  x-(3)  ...  x-(iV-A/  +  2)  I  ^^2) 


[  x(l)  •••  x{N-M)  x*(A/-|-l)  ...  x*(A^)  J 

and  b  is  the  desired  response  vector  given  by, 

■  x(A/-|-l)  ■ 

x{N) 

b=  •••  (4.3) 

x-(l) 

.  x’{N  -M)  . 

where  *  denotes  is  the  complex  conjugation  operator,  N  is  the  number  of  data  points 
and  M  is  the  desired  order  of  the  tap-weight  vector.  The  desired  response  vector  b  of 
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Eq  4.3  is  defined  as  the  next  data  point  of  the  cor  j^^ading  rows  of  .  The  right 
hand  side  (RHS)  of  Eq  4.1  can  be  written  as 

e  =  A^l.  (4.4) 


The  le2ist  squares  solution,  or  the  sum  of  th  error  squares,  is  minimized  only 
when  the  estimation  error  vector  is  orthogonal  to  i.  estimate  of  the  desired  response 
vector.  The  data  matrix  A^  contains  both  the  forw  .rd  (left  hand  partition  of  Eq  4.2) 
and  backward  (right  hand  partition  of  Eq  4.2)  cova;  ance  matrices,  which  doubles  the 
matrix  size.  Since  statistics  of  a  stationary  proces  are  assumed  equivalent  for  both 
directions,  the  solution  benefits  from  a  form  of  a-  ^raging.  The  uncorrelated  noise 
tends  to  average  to  a  lower  value,  while  the  correlated  signals  are  enhanced  in  the 
processing.  A  concise  explanation  of  linear  predi^  on  leading  to  the  pseudoinverse 
solution  is  presented  in  (Ref.  9]  and  is  summarized  here  for  clarity.  Only  the  forward 
covariance  data  matrix  is  used  in  the  developiucui.  fc,eneralization  to  the  forward- 
backward  data  matrix  of  Eq  4.2  is  straightforward. 

Linear  prediction  can  be  illustrated  using  a  sliding  window  [Ref.  10]; 


u;i  W2 


Xs  Xff^l 


xm 


X2  Xi 


The  dot  product  of  the  coefficient  array  w  with  :  .le  data  array  can  be  written  in 
matrix  form  as; 


Xm 

■  X2 

Xl 

UJ, 

a^M+1 

Xm+1 

■  X3 

X2 

, 

U>2 

= 

^M+7 

.  Xn^i 

.  .  .  . 

XS-M  . 

.  . 

.  . 

A  .  w  =  X 

data  matrix  tap  —  weight.'  data  predictions 


(4.5) 


where 


M  —  model  ord  r 
N  —  data  lengt  . 
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The  problem  is  uniquely  specified  when  the  data  matrix  in  Eq  4.5  is  square  and  of 
the  same  order  as  the  tap- weight  vector.  The  problem  is  overspecified  when  the  data 
matrix  is  larger  than  the  tap-weight  vector  and  the  number  of  solutions  is  given  by 
where  N  is  the  number  of  data  points  and  M  is  the  order  of  the  tap- weight 
vector.  The  overspecified  case  requires  least  squares  techniques.  Equation  4.5  can  be 
represented  as  a  discrete  difference  equation  and  the  presence  of  a  single  undamped 
sinusoidtil  frequency  component  is  then  described  by  a  conjugate  pole  pair  on  the 
unit  circle.  In  general,  2M  data  points  are  required  to  solve  for  M  sinusoids  as 
demonstrated  in  the  second  order  example, 


X2 

Xl 

in. 

Xz 

X3 

x-i 

u>2 

.  . 

The  pseudoinverse, 

ib  =  {A”A)-'A"b,  (4.6) 

accommodates  the  overspecified  set  of  linear  equations  without  compromising  the  so¬ 
lution  for  the  uniquely  determined  case.  A  deterministic  solution  exists  for  noiseless 
data  or  a  best  fit  in  the  least  squares  sense  exists  for  noisy  data.  The  only  difference 
in  Eq  4.1  and  Eq  4.6  for  this  application,  is  the  A^ A  matrix.  Equation  4.1  is  taken 
to  be  the  modified  covariance  matrix,  while  A^A  in  Eq  4.6  is  t2dcen  to  be  the  for¬ 
ward  covariance  matrix  only.  The  ability  of  the  pseudoinverse  to  handle  overspecified 
systems  makes  the  same  equation  applicable  in  both  cases. 

When  the  data  array  is  large  compared  to  the  order  of  the  process  (i.e.  the 
number  of  sinusoids),  forming  A^  A  from  all  the  data  produces  unmanageable  matri¬ 
ces.  Consider  64  data  szimples  from  which  a  model  order  of  five  is  selected.  Using  a 
forward-backward  matrix  data  arr<mgement,  A^^  will  be  of  order  (5  x  118)  and  A  will 
be  of  order  (118  x  5).  The  resulting  order  (5  x  5)  A^  A  matrix  is  quite  manageable  but 
forming  the  product  is  cumbersome.  A  larger  data  set  would  present  a  considerable 
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increase  in  computational  and  storage  requirements.  A  way  out  is  provided  by  the 
RLS  matrix  update  method.  This  method  minimizes  data  storage  and  computational 
requirements  needed  to  form  A  but  produces  identical  results. 

1.  RLS  UPDATE  EXTENSION 

The  RLS  update  method  starts  by  forming  the  A^ A  covariance  matrix 
from  the  minimum  required  number  of  data  points.  Each  subsequent  sample  then 
dynamically  improves  the  covariance  matrix  through  a  recursive  equation  update 
which  can  be  written  as, 

=  (yl"/!).-,  +  (A"  A):,  (4.7) 


where  (A^A)^  is  formed  by  summing  an  outer  product  of  the  last  row  of  the  forward 
partition  and  an  outer  product  of  the  last  row  of  the  backward  partition  of  /!„,  the 
data  matrix  for  n  points.  These  are  added  to  the  previous  covariance  matrix  to  form 
the  new  covariance  matrix.  This  recursion  permits  a  least  squares  solution  for  a  large 
data  set  without  requiring  matrix  multiplies  beyond  the  selected  order  of  the  process. 
The  method  is  illustrated,  for  a  second  order  case  with  a  forward-backward  data 


arrangement,  as  follows;  [Ref.  9,  page,  38] 
N  =  <:,M=i2: 


(A^Ah  = 


+  xj*  +  xj2 
+  ijxj  +  ijij 


*1X2  +  *2*3  +  *j*;  +  *s*i 
if  +  x’  +  x”  +  *;* 


AT  =  5  W  =  2  : 


{AffAh  = 


*2  +  *3  +  *2*  +  ®3* 
*1*2  +  12*3  +  *J*J  +  *J*I 


*1*2  +  *2*3 


,  +  xi*;  +  x;x; 

+*r+*;*  J^l 


*3*4  +*I* 


*3*4  +  *J*J 
*3  +  ^5 


/V  =  6  W  =  2  : 

=  f  *2  +  ^  +  *4  +  *5*  +  *1*  +  +  *7*3  +  *3*1  t  *I®S 

[  *1*2  +*2*3  +*3*4  +*;*J +*j*;  +  *i*;  xf +x*+x*+*”+x;*+x;* 


+ 


r  *f  +  *j*  *4xs  +  x*x;  ] 

[  X4XS  +  x*Xg  xf  +  xj*  J 

The  recursive  update  process  must  also  be  applied  to  the  RHS  of  Eq  4.1 


which  is  the  written  in  terms  of  6  in  Eq  4.4.  The  update  using  the  0  representation 
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can  be  written  as, 


=  ^n-1  +  (4.8) 

where  the  recursive  update  is  formed  by  summing  the  previous  value  with  an 
inner  product  of  a  matrix  consisting  of  the  last  column  of  the  forward  partition  and 
the  last  column  of  the  backward  partition  of  with  the  new  additions  to  desired 
response  vectors  of  6„  for  the  forward  and  reverse  partitions.  For  example, 


N  =  4: 

M  = 

2  : 

X3 

A»  = 

X3  xJ  Ij 

1  3ri 

^2  *4  . 

H 

t _ 

g  _  '  rjrs  +  1:3x4  +  ijx*  +  ijxj 


[  X1X3  +  X2X4  +  xjxj  +  xjxj  J 
N  =  5  M  =  2  : 

g  _  [  X3X3  +  X3I4  +  xjxj  +  xjxj  I  f  l4X5+xJxJ 

*  [  11X3  +  X,X4  +  xjxj  +  xjxj  j  [  X3I5  +  xjxj 

N  =  e  „  M  =  2 

g  _  f  XJX3  +X3X4  +X4X5  +xjxj  +i;xj  +xji;  I  XSI6  +  X*xJ  1 

■■  [  XIX3  +XaX4  +X3XS  +xjxj+xjxj+xjxj  J  X4X8+xJxJ  J  ‘ 

To  clarify  further,  the  matrix  for  the  (2  x  2)  example  above  where  =  5  is  as 
follows, 


N  =  5  :,  M  =  2  :  A» 


xj  xj  X4  .  x;  ij  x; 


XI  X3  X3  :  xj  i;  xj 


The  last  column  of  the  forward  partition  and  the  last  column  of  the  backward  partition 
along  with  the  new  points  to  be  predicted  from  the  b  vector  can  be  written  in  the 
update  matrix  equation  as. 


\  1  .  f  1  _  [  l4Xj+x;xJ  1 

I  ^3  J  I  J  ~  I  ^3XJ  +  ijxj  J 

which  is  the  matrix  addition  for  the  N  =  5  step  in  the  previous  example. 

The  recursive  update  for  the  covariance  matrix  and  the  desired  response 
vector  of  Eq  4.7  and  Eq  4.8,  allows  the  introduction  of  a  weighting  factor  that  can 
accommodate  a  nonstationary  (shifting  spectral  poles)  data  stream.  These  equations 
can  be  generalized  as; 

(A"A)„  =  X{A^A)n.i  +  {A^A)c,  (4.9) 


51 


and 


(4.10) 


where  A  is  an  exponential  weighting  factor  applied  to  past  data.  To  quote  from  [Ref. 
9,  page,  39],  “Typically,  A  is  less  than  unity,  thereby  aging  out  old  data,  hence  the 
expression.  Forgetting  Factor”  A  unity  forgetting  factor  would  simply  return  the 
least  squares  solution  of  Eq  4.1.  To  set  A,  simply  select  the  number  of  updates  (some 
time  interval  based  on  sampling  rate)  and  decide  the  level  of  suppression  required  at 
that  point.  For  example,  suppose  the  desired  suppression  is  10”®  after  60  samples. 
Then  A®®  =  10”®  and  A  =  0.8254. 

The  forgetting  factor  contributes  to  the  overall  numerical  stability  of  the 
algorithm,  when  it  is  employed  on  large  data  sets.  The  diagonal  elements  of  the 
covariance  matrix  are  sums  of  squares.  In  a  finite  arithmetic  implementation,  if  the 
data  are  not  zero  mean,  A^A  will  numerically  saturate,  since  normalization  of  the 
matrix  is  not  included  in  the  algorithm.  Even  if  the  data  are  of  zero  mean,  the 
numbers  along  the  diagonal  will  grow  steadily  if  A  is  unity.  A  forgetting  factor  less 
than  unity,  alleviates  the  problem  of  saturation  in  most  cases.  However,  any  data 
should  be  processed  with  the  mean  removed  to  minimize  the  likelihood  of  a  numericad 
problem. 


The  equation. 


(4.11) 


represents  the  algorithm  thus  far.  The  object  is  to  find  a  solution  for  u;„  of  the 
selected  order.  The  covariamce  matrix  (A^A)„  and  the  vector  6n  contain  a  description 
of  the  process  for  a  period  dictated  by  the  forgetting  faictor.  The  obvious  solution 
is  to  compute  [A^ A)~^.  This  is  computationally  inefficient  and  is  not  an  optimal 
approach.  Traditional  RLS  would  utilize  the  matrix  inversion  lemma  to  develop  a  set 
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of  recursive  equations  that  avoid  direct  computation  of  the  inverse,  yet  will  converge 
to  the  correct  result  after  some  number  of  iterations. 

2.  THE  MFBLP  SOLUTION 

At  this  point,  the  MFBLP  algorithm  described  in  [Ref.  6,  page,  349], 
can  be  used  with  minor  modifications  to  produce  high  resolution  dynamic  spectra 
of  the  process  from  the  {A^A)n  matrix.  Chapter  7  of  [Ref.  6]  provides  detailed 
derivations  and  discussions  of  least  squares  problem  and  the  reader  is  referred  there 
for  a  rigorous  treatment.  The  salient  points  are  summarized  here  as  they  apply  to 
the  implementation  employed.  The  subscript  n  has  been  inserted  in  the  equations 
presented,  to  emphasize  the  RL3  extension  that  will  be  utilized  in  the  final  combined 
algorithm.  The  subscript  is  to  be  suppressed  for  the  conventional  MFBLP  algorithm 
discussion  that  follows. 

Letting  Af  denote  {A^ A)~^ A^ ,  the  pseudoinverse  matrix,  Eq  4.11  can  be 
rewritten  as, 

Wn  =  A*bn  (4.12) 


where  all  elements  of  the  equation  have  been  previously  defined. 

A*  is  also  defined  in  terms  of  the  singular  value  decomposition  (SVD).  A 
statement,  without  derivation,  of  this  formulation  of  the  pseudoinverse  is  [Ref.  6], 


A*  =  X„ 


0 

0  0 


(4.13) 


where. 


E„  =  diag{<T^^  ,  ,  •  *  • ,  <^wn) 


and  the  cr„’s  are  the  singular  values  of  the  pseudoinverse  matrix.  Wn  is  the  rank  of 
the  matrix  and  if  A*  is  Hermitian,  the  singular  values  are  simply  the  absolute  value 
of  the  eigenvalues  of  A*.  In  general,  if  A*  were  an  order  Lx  M  matrix,  would  be 
&  M  X  M  matrix  of  columns  of  eigenvectors  of  {A^ A)n  and  Yj^  would  be  an  L  x  L 
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matrix  of  rows  of  eigenvectors  of  {AA^)n.  The  matrix  containing  would  be  order 


M  X  L. 


Substituting  Eq  4.13  into  Eq  4.12  yields, 


_  y  n  ” 


E-'  0 


Y”b 


which  can  be  partitioned  as, 


E:'  0  1  r  K, 


This  reduces  to 


=  x,„T:'y»K 


and  using  the  result  from  [Ref.  6,  Eq  7.79,  page,  333], 


then, 


U.,  =  XinS;‘xliA”i>„ 


which  in  terms  of  summations  becomes, 


(4.14) 


(4.15) 


(4.16) 


(4.17) 


(4.18) 


(4.19) 


Using  Eq  4.4  and  the  fact  that  is  an  eigenvector  of  the  deterministic  correlation 
matrix  {A^ A)n  with  associated  eigenvalue  A,„  =  the  critical  result  becomes, 


in  =  E 

1=1  ^‘*1 


(4.20) 


This  result  allows  computation  of  the  tap-weight  vector  of  an  AR  process  from  the 
eigenvalues,  eigenvectors  and  a  simply  computed  desired  response  vector.  The  sig¬ 
nificance  of  this  result  cannot  be  understated.  To  clarify  the  extensions  that  are 
employed  in  the  new  algorithm,  the  steps  used  in  the  conventional  MFBLP  method 


are  described.  These  are; 
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1.  Form  the  data  array  according  to  Eq  4.2  arrangement. 

2.  Use  the  SVD  on  to  compute  the  eigenvalues  and  eigenvectors  of  the  M  x  M 
deterministic  correlation  matrix  {A^ A). 

3.  Separate  the  eigenvalues  belonging  to  the  signal  and  those  belonging  to  the 
noise  subspaces.  The  signal  subspace  contains  the  K  dominant  eigenvectors. 

4.  Form  the  9  vector  according  to  Eq  4.4. 

5.  Use  0  and  the  signal  subspace  eigenvectors  and  eigenvalues  (Eq  4.20)  to  compute 
the  predictor  tap-weight  vector. 

6.  Form  the  (A/  -f  1)  x  1  prediction  error  filter  tap- weight  vector, 


a  = 


7.  Use  the  a  tap-weight  vector  to  compute  the  angulcir  frequency  of  the  of  the 
sinusoids  as  peaks  of  the  sample  spectrum  according  to, 


|a"5(u;)|^ 

where  s(u;)  is  the  (Af  -i-  1)  x  1  sinusoidal  signal  vector: 


1 

exp  —ju; 

_  exp  —jMu 


Tufts  and  Kumaresan  [Ref.  11]  have  experimentally  determined  the  opti¬ 
mal  order  for  this  algorithm  to  be  A/  =  37V/4.  In  its  present  form,  the  algorithm 
performs  well  as  a  frequency  estimation  tool,  but  to  quote  [Ref.  6,  page,  368],  “In  the 
conventional  FBLP  method,  5(u;)  represents  the  autoregressive  (AR)  power  spectrum 
of  the  process.  However,  this  is  not  so  in  the  modified  FBLP  method.”  Figure  4.1 
gives  an  example  spectra  computed  using  the  MFBLP  method  eis  described  above. 
The  64  point  test  data  set  was  obtained  from  [Ref.  12]  and  consists  of  two  closely 
spaced  equal  amplitude  sinusoids,  a  low  level  sinusoid  spaeed  away  from  the  dominant 
pair  and  colored  noise.  Spectral  truth  is  displayed  on  the  plot  in  dotted  lines  while  a 
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Figure  4.1:  Sample  MFBLP  spectra  of  the  Kay  test  data  set.  Eight  prin¬ 
ciple  eigenvalues/eigenvectors  and  48  coefficients  were  utilized  in  the  pro¬ 
cessing  (top  trace).  Conventional  periodogram  (lower  trace).  Spectral 
truth  (dotted  line). 


# 


.'>6 


standard  periodogram  of  the  data  set  is  also  included  using  a  dashed  line.  This  data 
set  will  be  used  for  tests  in  a  subsequent  section. 

Figure  4.1  clearly  demonstrates  the  basic  difficulty.  The  periodogram,  com¬ 
puted  using  a  rectangular  window,  has  problems  both  with  resolution  and  sidelobe 
interference.  Sidelobes  may  be  reduced  at  the  expense  of  resolution  but  little  else 
can  be  done  to  improve  the  spectral  estimate  using  standard  Fourier  processing.  The 
advantage  of  conventional  spectral  processing  is  that  relative  and  absolute  magnitude 
information  is  preserved. 

This  is  not  so  with  spectra  computed  using  the  MFBLP  technique.  The 
utilization  of  only  the  signal  subspace  eigenvalues  and  eigenvectors,  dictates  that 
spectra  computed  using  this  method  contain  only  the  principle  components  of  the 
time  series  (i.e.  the  correlation  matrix  does  not  reflect  the  complete  process).  It 
is  evident  from  Fig  4.1  (top  trace);  that  relative  and  absolute  magnitudes  are  not 
preserved.  The  figure  clearly  shows  the  accuracy  of  the  frequency  information  ioi 
the  dominant  components  and  the  lack  of  other  spectral  information.  Fortunately 
an  extension  exists  which  combines  the  benefits  of  periodogram  processing  with  the 
principle  component  enhancement  of  the  MFBLP  technique. 

3.  SPECTRA  USING  LINEAR  PREDICTION 

Instead  of  relying  on  the  spectra  from  the  tap-weights  to  characterize  the 
process,  the  weights  are  used  to  extend  the  data  via  linear  prediction.  Conventional 
periodogram  analysis  is  then  used  on  the  extended  data  set.  This  modification  to  the 
conventional  MFBLP  technique  is  investigated  in  detail  in  [Ref.  8]  with  impressive 
results.  The  reader  is  referred  to  this  thesis  for  a  comparison  of  the  MFBLP  method 
with  many  other  modern  spectral  estimation  techniques.  Its  performcince  is  noticeably 
superior  in  many  respects.  Figure  4.2  demonstrates  the  performance  of  the  technique 
on  the  Kay  data  set.  The  64  point  data  set  is  extended  in  the  forward  and  backward 
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directions  by  96  points,  yielding  a  new  time  series  of  256  samples.  The  important 
parameters  used  for  the  extension  were,  48  tap-weights  computed  using  the  first 
eight  principle  eigenvalues  and  eigenvectors  of  the  covariance  matrix  decomposition. 
Rectangular  windowed  and  Hamming  windowed  periodograms  of  the  original  64  point 
time  series,  zero  padded  to  512  points,  have  been  overlayed  for  comparison,  along  with 
spectral  truth.  Note,  to  preserve  clarity  in  the  plot  a  windowed  version  of  the  extended 
periodogram  was  not  included.  The  double  peak  observed  for  the  low  frequency  low 
level  signal  has  a  much  improved  character  when  a  Hamming  window  is  applied  to 
the  extended  series,  as  can  be  observed  in  subsequent  plots. 

Success  of  the  technique  can  be  attributed  to  the  extension  of  the  principle 
components  in  the  data.  The  added  points  increase  the  real  resolution  of  the  FFT, 
but  more  importantly  the  sidelobe  structure  is  much  improved  and  low  level  signals 
present  in  the  data  become  visible  eventhough  information  about  them  is  not  present 
in  the  tap-weights.  Standard  windowing  techniques  can  also  be  apphed  to  improve 
sidelobe  structure  in  the  spectral  estimates  if  desired.  Both  relative  and  absolute 
magnitude  relationships  are  preserved  relatively  well  with  this  method.  In  many 
applications  this  is  an  important  consideration. 

At  this  point,  operation  of  the  combined  algorithm  should  be  reasonably 
evident.  Briefly  the  combined  method  will  operate  as  follows.  The  modified  co- 
variance  matrix  {A^ A)n  is  formed  from  the  data  using  the  forward-backward  data 
arrangement.  This  matrix  is  updated  using  the  RLS  update  technique  discussed 
earlier.  At  any  update  a  spectrum  can  be  produced  using  the  extended  version  of 
the  MFBLP  technique  cis  described  earlier  in  this  section.  An  important  difference 
between  the  combined  technique  and  the  MFBLP  is  that  an  SVD  is  performed  on 
the  An  or  A^  matrix  for  the  MFBLP  technique.  In  the  combined  technique,  this 
matrix  is  not  available.  An  eigenvalue  decomposition  is  performed  on  the  covariance 
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Figure  4.2:  Sample  spectra  of  the  MFULP  technique  with  the  linear  pre¬ 
dictive  extension,  applied  to  the  Kay  test  data  set. 


f 
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matrix,  {A^A)n-  The  MATLAB  SVD  algorithm  can  provide  the  eigenvalue  decom¬ 
position  of  the  covariance  matrix  so  it  is  utilized  for  its  numerical  robustness.  Any 
eigenvalue  decomposition  method  could  be  used  without  altering  the  result.  The 
eigenvalue  decompositions  are  the  only  computationally  expensive  aspect  of  the  total 
algorithm.  Only  the  principle  eigenvalues  and  eigenvectors  are  required  for  operation 
of  the  algorithm,  so  any  efficient  method  for  estimating  the  partial  decomposition  of 
the  covariance  matrix  to  get  the  principle  components  could  be  utilized  to  improve  ef¬ 
ficiency.  Algorithms  for  this  purpose  have  been  discussed  in  the  recent  literature.  The 
term  SVD  used  from  this  point  on  will  refer  to  the  MATLAB  SVD  algorithm.  The 
MFBLP  method  described  in  [Ref.  6]  is  the  basis  for  two  slightly  different  extensions 
that  yield  improved  performance.  For  lack  of  better  temunology  but  to  differenti¬ 
ate  between  the  three  techniques  discussed,  the  combined  algorithm  has  been  labeled 
modified  recursive  least  squares  (MRLS).  To  emphasize  the  changes  applied  to  the  con¬ 
ventional  MFBLP  technique,  it  will  be  termed  extended  modified  forward-backward 
linear  prediction  (XMFBLP). 

B.  OUTLINE  OF  MRLS  SPECTRAL  ESTIMATION 

It  follows  from  previous  discussions  that,  with  a  forgetting  factor  of  unity  in  the 
RLS  update  (Eq  4.7  and  Eq  4.8)  and  a  fixed  data  length,  MRLS  and  XMFBLP  will 
produce  identical  spectral  results.  Although  it  is  normal  for  the  XMFBLP  to  utilize 
the  whole  data  set  at  once,  it  is  perfectly  acceptable  to  form  the  data  matrix  on  any 
contiguous  subset  of  data  arranged  to  yield  the  desired  processing  order.  One  could 
conceive  an  algorithm  that  moved  through  a  large  data  set  point  by  point,  turning 
out  the  first  point  in  the  data  matrix  as  it  accepts  a  new  point  in  the  last  position. 
This  approach  would  be  akin  to  segmented  periodogram  processing  with  overlap  and 
the  spectra  produced  would  be  local  high  resolution  snapshots  of  the  process  for  each 
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new  data  point.  Since  the  output  for  a  data  point  would  include  a  covariance  matrix 
computed  with  only  the  present  data  segment,  long  term  trends  or  low  level  periodic 
signals  might  be  missed  without  some  conventional  spectral  averaging.  An  advantage 
however,  is  that  no  ambiguity  would  exist  with  respect  to  where  to  apply  the  linear 
predictive  data  extension  of  the  XMFBLP  algorithm.  It  would  automatically  occur 
at  the  ends  of  the  data  subset  being  analyzed. 

By  contrast,  the  RLS  update  method  provides  for  successive  improvement  of 
the  covariance  matrix  as  each  new  data  point  is  added.  The  period  over  which  the 
matrix  is  valid  is  dependant  on  the  value  of  the  forgetting  factor  over  all  the  data 
from  the  start  of  processing.  Because  the  covariance  matrix  can  be  made  to  reflect 
trends  over  a  much  longer  period  than  the  order  of  the  matrix,  the  points  at  which  the 
linear  predictive  data  extension  should  be  applied,  to  generate  a  spectral  estimate, 
are  somewhat  ambiguous.  The  successive  improvement  of  the  covariance  matrix  af¬ 
forded  by  the  RLS  update  method,  should  provide  superior  performance  for  low  level 
periodic  signals  in  a  large  data  set  without  the  need  for  averaging  of  the  individual 
periodograms,  although  this  remains  ah  additional  processing  option.  For  dynamic 
spectral  processing  of  nonstationary  data,  it  is  the  tracking  of  the  shifting  poles  of  the 
process  that  is  the  desired  measurement  and  the  use  of  the  MRLS  algorithm  provides 
an  excellent  opportunity  to  provide  reasonable  accuracy  in  this  respect. 

The  purpose  of  the  data  extension,  as  described  earlier,  is  not  to  accurately 
predict  the  next  data  point  in  the  sequence.  The  next  point  can  either  be  measured, 
or  is  already  available.  The  intent  is  to  sharpen  the  spectral  frequency  and  magnitude 
measurements  of  the  principle  components  at  work  in  the  data  at  a  particular  time 
cind  present  them  in  a  high  resolution  display,  so  changes  over  time  can  be  observed. 
Logically  then,  for  the  MRLS  method,  the  data  extension  is  apphed  to  a  data  subset 
that  consists  of  the  amount  of  data  that  would  be  utihzed  by  the  same  order  XMFBLP 
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method  if  it  were  applied  at  the  present  update  sample.  This  is  not  necessarily 
optimum,  but  no  inaccuracy  is  induced. 

Steps  in  the  MRLS  algorithm  may  be  summarized  ais  follows: 

1.  Form  the  modified  covariance  matrix,  {A^ A)i,  and  the  vector  (Eq  4.4)  from 
the  minimum  data  required  to  meet  the  desired  model  order  M. 

2.  Select  the  value  of  the  forgetting  factor,  A,  and  update  Eq  4.7  and  Eq  4.8 
recursively,  for  each  successive  data  sample. 

3.  Compute  a  spectra  at  any  point  “n”  or  desired  interval  by  performing  an  eigen¬ 
value/eigenvector  decomposition  to  the  covariance  matrix,  {A^ A)n. 

4.  Select  the  signal  subspace  eigenvalues  and  eigenvectors.  The  inclusion  of  some 
eigenvectors  and  eigenvalues  from  the  noise  subspace  will  not  greatly  affect  the 
results,  so  the  number  selected  can  remain  constant  throughout  a  processing 
run. 

5.  Apply  Eq  4.20  to  compute  the  tap- weight  vector. 

6.  Use  the  tap-weight  vector  as  a  linear  predictive  filter  to  extend  the  local  time 
series  forward  and  backward  as  described  previously  in  this  section. 

7.  Apply  a  conventional  periodogram  (windowing  optional)  to  the  extended  data 
subset,  to  compute  the  local  spectrum  for  the  particular  period. 

8.  Save  the  individual  spectra  to  produce  a  dynamic  spectral  display  of  the  process. 

C.  MRLS  TESTING  AND  PERFORMANCE 

A  sample  test  data  set  w<is  compiled  by  concatenating  severiil  64  point  specific 
test  cases  into  a  single  time  series.  Each  test  case  in  the  time  series,  was  separated 
from  the  next  test  case  by  64  points  of  white  gaussian  noise.  The  forgetting  factor  was 
set  to  ensure  that  the  memory  window  of  the  covariance  matrix  was  approximately 
64  points.  Real  data  would  more  than  likely  have  smooth  transitions  rather  than  the 
step  transitions  of  this  test  data  series,  and  thus,  the  test  series  provided  a  formidable 
test  for  the  algorithm.  Table  4.1  contains  a  list  of  the  characteristics  of  each  test  data 
subset  in  the  total  time  series. 
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TABLE  4.1:  A  TABLE  OF  SIGNALS  COMPRISING  THE  MRLS  TEST 
DATA  SET. 


Description 

#  Points 

C  haracteristics 

Test 

Gaussian  noise 

64 

Unity  variance 

Algorithm  initialization 

Gaussian  noise 

64 

Unity  variance 

Noise  alone 
performance 

Kay  data  set 

64 

3-sinusoida,  colored  noise 

fi  =  6.4Hz,  /j  =  12.8Hz,  fi  =  13.44Hz 

Resolution 

performance 

Gaussian  Nois« 

64 

Unity  variance 

Case  separator 

Multiple  signals 

64 

3-sinusoids  SNR=3dB 

fi  =  5.0Hz,  h  =  lO.OHz,  /3  =  IS.OHz 

Equal  amplitude 
performance 

Gaussian  Noise 

64 

Unity  variance 

Case  separator 

Multiple  signals 

64 

3-sinusoids  SNR=-3dB 

/i  =  S.OHz,  f2  =  lO.OHz.  h  =  IS.OHx 

SNR 

performance 

Gaussian  Noise 

64 

Unity  variance 

Case  separator 

Gaussian  window 

64 

l-binusoid  SNR=:5dB 
/  =  IS.OHz 

Variable  amplitude 
performance 

Gaussian  Noise 

64 

Unity  variance 

Case  separator 

Variables  in  the  MRLS  processing  were  set  cis  follows; 

1.  The  order  selected  was  47,  which  is  approximately  3 A/4  where  N  is  64.  The 
cases  of  interest  are  contained  in  64  point  blocks,  and  this  is  the  experimentally 
determined  optimum  for  the  MFBLP  as  described  earlier. 

2.  The  forgetting  factor  was  set  to  be  A  =  0.9  which  corresponds  to  a  suppression 
of  approximately  0.001  after  64  updates. 

3.  The  data  is  extended  forward  and  backward  96  samples  in  each  direction  for 
each  spectral  computation.  The  new  data  length  of  256  points  is  zero  padded 
to  512  points  for  the  periodogram. 

4.  A  Hamming  window  was  applied  before  periodogram  processing. 

5.  The  principle  eigenvalues  and  eigenvectors  were  taken  to  be  the  first  six  from 
the  eigenvailue  decomposition. 


The  algorithm  was  initialized  on  the  first  64  points  of  gaussian  noise  data  and 
was  allowed  to  run  on  the  full  640  point  data  set.  Outputs  were  taken  before  each  case 
separation  noise  sequence.  For  comparison,  these  outputs  are  overlayed  with  conven¬ 
tional  periodograms  of  the  64  point  test  cases  zero  padded  to  512  points,  MFBLP 
plots  of  the  tap- weight  spectra  and  spectra  computed  using  the  XMFBLP  technique. 
Spectral  truth  is  also  included  on  each  plot.  The  point  in  the  time  series  at  which 
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Point  #  128,  max  4.96934 


Figure  4.3:  Noise  alone  performance  of  the  MRLS  algorithm. 

the  spectra  are  taken  and  the  maximum  numerical  value  in  the  covariance  matrix  at 
that  time  are  displayed  at  the  top  of  each  plot. 

Figure  4.3  shows  the  noise  only  performance  of  this  algorithm.  All  noise  pro¬ 
cesses  contain  some  structure  especially  when  the  sequences  are  short  duration.  The 
figure  shows  that  the  absolute  level  of  the  spectrum  is  depressed  in  comparison  to 
the  conventional  pcriodogram  of  the  test  case,  which  is  expected.  The  structure  is 
somewhat  enhanced  (i.e.  peaks  are  slightly  more  pronounced)  but  there  is  close  cor- 
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respondence  between  the  conventional  periodogram  and  MRXS  results.  The  noise 
suppression  is  a  result  of  the  additional  resolution  and  related  distribution  of  the 
original  noise  content  over  more  frequency  bins.  T.*'e  structure  is  no  worse  than  the 
conventionaJ  periodogram  result  and  provides  an  improvement  in  absolute  terms.  The 
whiteness  of  the  noise  would  be  more  evident  if  successive  reahzations  were  averaged. 

Figure  4.4  demonstrates  the  performance  on  the  Kay  data  test  series  described 
earlier.  Both  the  MRLS  spectrum  and  the  conventional  periodogram  have  Hamming 
windows  applied.  The  performance,  in  this  case,  is  superior  in  most  respects  to  the 
conventional  periodogram  result.  The  two  closely  spaced  signals  are  easily  resolved 
and  at  the  correct  spectral  levels  considering  the  Hamming  window  affects  on  the 
sinusoids.  The  low  level  signal  is  visible  but  is  depressed  from  its  actual  value.  The 
tap- weight  spectrum  or  conventional  MFBLP  result  shows  no  enhancement  of  the 
low  level  signal  thus  it  suffers  the  same  degradation  as  the  noise.  Also  note  that 
because  there  is  no  enhancement  of  this  component,  it  does  not  enjoy  the  same  peak 
resolution  as  the  higher  amplitude  signals  but  it  is  still  better  than  the  conventional 
periodogram.  The  narrowed  sidelobe  structure  over  that  of  the  conventional  peri¬ 
odogram  is  also  a  useful  feature  of  the  method.  The  signal  peaks  do  not  reach  the 
magnitude  levels  shown  by  the  truth  hues  because  of  the  broadening  of  the  bin  main 
lobes  and  subsequent  spreading  of  the  sinusoidal  energy  due  to  the  application  of  the 
Hamming  window.  Other  processing  with  a  rectangular  window  has  shown  the  peak 
levels  of  the  high  amplitude  signals  match  the  truth  levels.  The  Hamming  windowed 
conventional  periodogram  appears  to  better  reflect  the  shape  of  the  colored  noise  part 
of  the  spectrum  in  Fig  4.4.  This  is  probably  because  the  bin  width  is  large  enough  to 
mask  the  noise  sub-structure  in  this  short  realization.  The  MRLS  algorithm  enhances 
various  principle  peaks  in  the  colored  noise  region  which  gives  an  indication  of  what 
might  actually  be  occurring  in  that  part  of  the  spectrum  for  this  realization. 
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Point  #  320,  max  17.6022 


Figure  4.5:  MRLS  performance  with  3.0  dB  SNR  and  multiple  sinusoids. 

Figure  4.5  and  Fig  4.6  give  an  indication  of  the  SNR  performance  of  the  al¬ 
gorithm  with  multiple  signals  of  equal  amplitude.  Although  most  modern  methods 
require  greater  than  10  dD  of  SNR  ratio,  this  method  performs  very  well  at  3  dB 
SNR  and  still  shows  good  response  at  a  -3  dB  SNR.  There  is  some  frequency  bias, 
but  the  figures  demonstrate  the  performance  increase  over  conventional  periodogram 
processing. 

Figu  res  4.7,  Fig  4.8  and  Fig  4.9  give  an  indication  of  the  dynamic  performance 
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Figure  4.6:  MRLS  performance  with  -3.0  dB  SNR  and  multiple  sinusoids. 
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of  the  aigorithm  on  a  sinusoidal  packet.  A  15  Hz  sinusoid  is  multiplied  by  a  Hamming 
window  yielding  a  5  dB  SNR  at  the  center  of  the  packet.  Figure  4.7  is  16  points  prior 
to  the  entire  64  point  signal  updating  the  covauiance  matrix.  Figure  4.8  shows  the 
spectral  result  after  the  matrix  is  updated  with  all  64  points  and  Fig  4.9  shows  the 
effect  of  an  additional  16  updates  of  gaussian  noise.  Curiously,  MRLS  appears  to 
perform  better  at  a  plus  and  minus  16  points  than  it  does  with  the  maximum  number 
of  process  samples  having  updated  the  covariance  matrix.  The  MFBLP  tap-weight 
spectra  shows  very  little  enhancement  for  this  case,  thus  the  entire  MRLS  spectrum 
is  depressed.  Using  more  eigenvalues  and  eigenvectors  in  the  generation  of  the  tap- 
weights  would  likely  change  this  result. 

D.  MRLS  SUMMARY 

The  MRLS  algorithm  was  developed  to  provide  a  reasonable  compromise  of 
signal  processing  parameters  to  track  spectral  information  of  an  unknown  process. 
The  algorithm  provides  reasonable  SNR  performance,  excellent  frequency  resolution 
capability,  a  point  by  point  process  frequency  tracking  capability,  and  numerical 
stability.  The  data  set  for  which  it  is  intended  is  comprised  of  six  hour  data  sections 
yielding  approximately  11000  sample  points  spaced  1.9375  seconds  apart.  In  order  to 
look  at  the  lower  frequencies  the  data  set  must  be  filtered  and  decimated.  K  the  data 
set  is  decimated  by  10  then  the  number  of  points  available  for  processing  drops  to 
1100.  This  is  not  a  lot  of  points  to  determine  spectral  dynamics,  thus  the  necessity 
to  implement  a  higher  resolution  nonstationary  spectral  estimation  scheme.  Results 
of  processing  the  tomography  experimental  data  set  are  presented  in  Chapter  V. 

The  advantages  of  the  MRLS  processing  scheme  eire  summarized  as  follows: 

1.  Simple  update  procedure. 

2.  No  Matrix  Inverse  required. 

3.  Stable  numerical  techniques  utilized. 
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Figure  4.7:  MRLS  performance  on  a  Hamming  windowed  sinusoidal  packet 
burst  with  5dB  of  SNR  at  the  packet  peak  value,  -16  points. 
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Point  #  576,  max  1.31427 


Figure  4.8:  MRLS  performance  on  a  Hamming  windowed  sinusoidal  packet 
burst  with  5dB  of  SNR  at  the  packet  peak  value,  all  test  points. 


Point  #  592,  max  8.0888 


Figure  4.9:  MRLS  performance  on  a  Hamming  windowed  sinusoidal  packet 
burst  with  5dB  of  SNR  at  the  packet  peak  value,  4-16  points. 
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4.  Excellent  noise  performance  for  a  modem  technique. 

5.  Handles  nonstationary  data  gracefully. 

6.  Excellent  resolution  performance. 

7.  Constant  iterative  improvement  of  the  covariance  matrix. 

8.  Low  memory  and  reasonable  computational  requirements  allow  it  to  be  imple¬ 
mented  with  good  performance  on  a  personal  computer  (PC),  eventhough  the 
data  set  may  be  quite  large. 

The  one  major  disadvantage,  is  that  the  present  implementation  requires  a  full 
covariance  matrix  eigenvalue  decomposition  for  each  spectral  output. 
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V.  RESULTS  AND  RECOMMENDATIONS 


The  scope  of  this  thesis  work  did  not  include  the  physical  interpretation  of 
data  processed.  The  intent  here,  was  to  develop  algorithms  that  w’ill  be  effective  in 
providing  spectral  information  in  the  surface  and  internal  wave  frequency  regions  of 
the  Hadamard  transform  matched  filter  output  of  the  MBTE  data  set.  This  section 
will  discuss,  in  general  terms,  the  outputs  of  the  LMS  tracker  for  arrivals  A,  B  and 
C  of  station  J  contained  in  Appendix  C  and  Appendix  D.  Conventional  periodogram 
analysis  is  used  to  verify  the  presence  of  the  surface  wave  component  in  the  error 
output  of  the  LMS  tracker.  It  would  interesting  to  investigate  the  dynamics  of  the 
surface  wave  spectra  using  the  MRLS  technique,  however,  internal  wave  spectral 
dynamics  are  of  more  interest,  so  results  of  processing  with  internal  waves  in  mind 
are  presented  instead. 

A.  STATION  J  ARRIVAL  TRACKING 

The  three  arrivals  traicked  in  the  station  J  data  have  been  defined  earlier  in  this 
thesis.  Plots  of  adl  the  outputs  of  the  LMS  tracking  routine  are  contained  in  Appendix 
C.  These  figures  are  arranged  in  groups  of  three.  Results  for  each  of  the  seven  outputs 
of  the  LMS  arrival  tracking  filters  are  presented.  This  provides  for  an  easy  comparison 
between  figures.  The  predicted  output  tracks  are  overlayed  on  the  correlograms  in 
Appendix  D.  Appendix  D  is  basically  a  set  of  truth  plots  to  determine  the  points  at 
w’hich  the  tracking  may  be  questionable  due  to  a  lack  of  signal.  No  such  validation  was 
available  in  earlier  processing.  Ninety-six  coefficients  were  used  for  processing  of  all 
three  arrivals  with  the  LMS  algorithm.  The  /i  parameters  0.003,  0.0015  an  0.009  were 
used  for  arrivals  A,  B  and  C  respectively.  The  small  differences  in  the  adaptation 
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Figure  5.1:  Track  comparison  for  Arrivals  A,  and  C  of  station  J. 

parameter  maintain  similar  variances  on  the  LMS  tracks  by  compensating  for  the 
differences  in  the  mean  of  the  arrival  levels. 

Figure  5.1  compares  the  LMS  tracks  for  the  station  J  six  hour  data  block  pro¬ 
cessed.  The  strongest  continuous  arrival  obvious  from  the  correlograms  is  arrival  B. 
The  performance  of  the  tracker  should  be  best  for  this  arrival  and  indeed,  the  cor- 
rclograms  verify  this  point.  In  Fig  5.1  the  arrival  A  and  C  tracks  show  some  large 
transitions.  Arrival  B  shows  only  one  smaller  transition.  The  correlograms  indicate 


that  this  variability  is  a  result  of  signals  that  fade  out  or  are  nonexistent.  When  a 
peak  is  lost,  the  algorithm  searches  for  another  peak  to  lock  on  to. 

Arrivals  A  and  C  have  unstable  sections  at  the  track.  The  stable  section  for 
arrival  A  begins  at  21.5  hrs  and  carries  to  the  end  of  the  track.  The  stable  section 
for  C  is  between  approximately  20.25  and  22  hrs.  Note,  the  time  scale  is  in  decimal 
hours.  All  the  plots  of  Appendix  C  show  hours  into  the  experiment  from  a  decimal 
start  time  thus  25  hrs  is  01:00  hrs  the  next  day  in  real  time.  Arrival  B  is  steady 
over  the  whole  tracking  period.  The  transition  in  arrival  B  appears  on  correlogram 
Fig  D.8.  There  are  some  complete  data  dropouts  on  this  correlogram  w’hich  may 
have  triggered  the  otherwise  stable  track  to  shift  peaks.  The  transition  might  also  be 
result  of  the  physical  processes  at  work  in  the  data,  as  it  stays  within  the  arrival  field 
and  it  is  quite  stable  on  either  side  of  the  transition. 

Looking  closely  at  the  arrival  B  character  on  the  correlograms,  one  can  discern 
a  periodic  amplitude  fluctuation  in  the  arrival.  Alternate  areas  of  small  dark  and 
light  packets  can  be  observed  to  occur  with  a  period  on  the  order  of  seconds.  These 
amplitude  swings  are  quite  large  and  are  hkely  caused  by  the  multipath  interference 
effects  described  earlier.  It  is  unfortunate  that  the  A  and  C  arrivals  were  not  more 
stable,  as  a  direct  visual  comparison  might  yield  obvious  points  of  general  similarity 
between  the  tracks  of  Fig  5.1.  There  is  only  a  very  short  region  where  all  three  tracks 
are  operating  on  good  SNR  data.  This  occurs  between  21  hrs  and  22  hrs  and  is 
not  large  enough  to  observe  the  trends  in  presentations  such  as  Fig  5.1.  However, 
the  scale  of  the  correlograms  does  show  similar  track  characteristic  for  the  B  and 
C  arrivals  around  21:50  hrs.  The  arrival  A  track  does  not  show  the  same  dynamics 
in  this  region  but  its  arrival  path  might  not  be  sampling  the  same  processes  as  the 
other  two  arrival  paths.  The  tracks  of  Fig  5.1  and  the  correlograms  produce  much 
information  useful  for  the  data  interpretation  but  other  LMS  outputs  also  provide 
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very  useful  information  as  well. 

Besides  arrival  time,  the  LMS  tracker  has  six  other  outputs  all  of  which  provide 
information  for  data  interpretation  and  performance  monitoring.  Gross  features  of 
these  outputs  provide  visual  indications  of  the  process  and  statistical  processing  of 
these  outputs  yields  much  additional  information.  For  instance,  the  arrival  time  error 
output  is  a  zero  mean  output  and  LMS  processing  attempts  to  keep  this  signal  as 
spectrally  white  as  possible.  The  variance  of  this  output  can  be  used  to  compare  di¬ 
rectly  with  modeled  variances  for  the  arrivals.  Also,  periods  when  this  signal  deviates 
markedly  from  zero  mean  indicate  regions  where  the  track  should  be  validated.  For 
example  Fig  C.9  shows  one  such  event  around  23  hrs  for  arrival  C.  A  check  of  the 
correlograms  around  this  time  show  that  tracker  has  lost  the  track  and  is  searching 
for  a  new  peak  to  lock  on  to.  Correlograms  show  that  it  does  not  succeed  in  regaining 
a  valid  track  after  this  point. 

The  amplitude  plots  of  Appendix  C,  which  include  predicted  amplitude,  mea¬ 
sured  amplitude  and  amplitude  error  for  all  three  arrivals,  show  some  slower  trends 
in  the  LMS  filtered  output,  but  more  importantly  the  variance  of  the  error  signals  are 
quite  large  for  all  three  tracks.  This  is  another  indication  of  the  interference  of  the 
closely  spaced  arrivals  in  the  data.  Arrival  B  has  the  most  stable  arrival  time  ♦rack 
with  the  lowest  arrival  time  error  variance.  This  can  be  observed  by  comparing  the 
arrival  time  error  signals  in  Appendix  C.  A  similar  comparison  of  the  amplitude  error 
signal  shows  that  arrival  B  has  the  largest  amplitude  variance  of  the  three.  Presum¬ 
ably,  if  arrival  B  is  a  sequence  of  single  fully  resolved  arrivals  then  it  would  have  the 
lowest  amplitude  variance  as  well,  since  the  amplitude  variations  are  more  accurately 
predicted  by  the  LMS  amplitude  tracker. 

The  phase  outputs  from  the  LMS  filtering  process  could  also  provide  useful 
information  for  the  interpreting  the  data.  The  phase  component  in  this  data  is  gen- 


erated  from  quadrature  sampling  of  the  original  time  series  before  matched  filtering. 
Both  components  of  the  quadrature  sampling  after  matched  filtering  are  used  to  com¬ 
pute  amplitude  and  phase.  The  amplitude  values  are  used  in  the  processing  described 
in  this  thesis.  The  relationship  of  the  computed  phase  to  the  phase  of  the  physical 
acoustic  signal  is  unclear.  The  lack  of  absolute  travel  time  information  contributes  to 
the  ambiguity  as  phase  unwrapping  without  this  information  might  not  be  possible. 
The  phase  tracks  output  do  show  trackable  behavior.  Despite  the  high  variance  a 
distinct  average  trend  is  evident  in  the  three  phase  plots.  Future  work  will  attempt 
to  exploit  the  phase  information  for  improved  tracking. 

The  seven  outputs  of  the  LMS  tracking  algorithm,  examples  of  which  are  in- 
clu  ied  in  Appendix  C  for  all  three  station  J  arrivals,  provide  a  rich  analysis  set  to 
which  many  forms  of  post  processing  can  be  applied  to  extract  information  useful  in 
the  physical  interpretation  of  the  experimental  data  set. 

B.  STATION  J  SPECTRAL  PROCESSING 

The  spectral  processing  was  performed  in  two  spectral  regions.  The  surface 
wave  region  and  the  internal  wave  region.  This  processing  uses  two  of  the  outputs 
from  the  LMS  tracks  described  in  the  previous  section.  For  spectral  estimation  the 
predicted  arrival  time  output  of  the  filter  is  used  for  the  internal  wave  region  and  the 
arrival  time  error  signal  is  used  for  the  surface  wave  region.  The  desirable  output 
in  this  section  is  a  spectral  output  that  can  show  the  dynamics  of  the  underlying 
processes.  The  primary  objective  for  the  surface  wave  region  in  this  thesis  was  to 
verify  its  presence  in  the  LMS  tracker  output.  This  was  done  using  conventional 
periodogram  techniques  on  the  arrival  time  error  signal.  The  results  are  summarized 
in  Fig  5.2  for  all  three  arrivals  and  are  comparable  with  the  wave  buoy  results  from 
Fig  5.3[Ref.  1].  The  full  10800  data  points  were  processed  using  non-overlapped  128 
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point  Hamming  windowed  peiiodogram  analysis.  The  resulting  periodograms  were 
averaged  to  produce  Fig  5.2. 

Figure  5.2  shows  that  arrival  B  has  the  closest  agreement  with  the  wave  buoy 
data  of  Fig  5.3.  Arrival  B  shows  a  double  peak.  It  is  possible  that  the  broadness 
of  the  wave  buoy  peak  is  due  to  the  presence  of  a  second,  lower  amplitude,  peak. 
The  resolution  of  the  arrival  B  processing  is  twice  that  of  the  wave  buoy  data  and 
the  averaging  is  three  time  as  long,  thus  it  is  better  suited  to  reveal  such  details. 
The  A  and  C  arrivals  of  Fig  5.2  show  slightly  different  character  than  arrival  B.  No 
attempt  has  been  made  to  separate  the  valid  track  sections  of  arrivals  A  and  C,  thus 
there  is  contamination  from  the  large  transitions  noted  earlier  in  the  tracks.  This  is 
likely  the  source  of  the  high  frequency  peaks  present  in  the  A  and  C  tracks  but  not 
present  in  the  cleaner  arrival  B  track.  This  contamination  would  also  explain  the 
higher  overcill  levels  of  the  A  and  C  spectra.  Note  despite  the  contamination  from  the 
track  transitions,  arrival  C  shows  a  reasonable  surface  wave  peak.  Intriguingly,  the 
peak  is  absent  in  the  arrival  A  spectrum.  Further  analysis  using  the  MRLS  technique, 
developed  in  Chapter  FV,  would  aid  in  explaining  the  character  of  these  spectra  and 
would  produce  dynamic  spectral  results  over  the  entire  time  period.  However,  having 
verified  the  presence  of  the  surface  wave  component,  this  is  left  in  favor  of  producing 
the  dynamic  spectra  in  the  internal  wave  frequency  region  using  the  MRLS  spectral 
estimation  technique. 

Figure  5.4,  Fig  5.5  and  Fig  5.6  show  the  results  of  the  dynamic  spectral  process¬ 
ing  in  3-D  waterfall  displays.  The  LMS  predicted  arrival  time  tracks  were  low  pass 
filtered  with  an  eighth  order  Chebychev  filter.  The  passband  had  a  corner  frequency 
of  0.02066  Hz  and  the  track  was  decimated  by  a  factor  of  10  after  the  low  pass  fil¬ 
tering.  These  decimated  tracks,  consisting  of  approximately  1100  data  points  each, 
were  input  to  the  MRLS  algorithm.  Spectra  were  plotted  at  each  16  point  update 
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Figure  5.2:  Spectra  at  surface  wave  frequencies  from  the  LMS  arrival  time 
error  tracks  for  arrivals  A,  B  and  C  of  station  J. 
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Sea  Surface  Spectrum 
NDBC  Buoy  14DEC88  2100  PST 


Frequency  (Hz) 


Significant  Wave  Height  3.54  m 
Average  Period  9.11  sec 
Dominant  Period  12.50  sec 
Dominant  Direction  314  N 


Figure  5.3:  Surface  wave  power  spectrum  in  Monterey  Bay  from  the  wave 
buoy  southwest  of  Santa  Cruz,  14  Dec  88.  This  spectrum  is  computed 
from  two  hours  of  data  ending  at  2100  hrs  PST. 
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of  the  order  47  covariance  matrix  used  in  the  MRLS  processing.  Other  M!;LS  pa¬ 
rameters  included  a  forgetting  factor  of  0.9  with  six  principle  eigenvalues  employed, 
a  sub-series  extension  of  96  points  in  both  the  forward  and  backward  directions  and 
a  Hamming  window  applied  in  the  periodogram  computation.  The  spectra  were  con¬ 
verted  to  decibels  for  display  in  the  3-D  waterfalls.  Each  spectr?’  trace  shown  in  the 
dynamic  plots  represents  approximately  5  minutes  of  data.  This  resolution  could  be 
increa.sed  by  a  factor  16.  For  illustration  purposes  the  5  minute  resolution  is  sufficient. 
Additionally,  the  decimation  yields  nine  other  realizations  that  could  be  averaged  if 
required. 

The  results  in  these  figures  are  very  interesting.  There  is  no  way  to  verify  the 
presence  of  internal  waves  in  this  data  set  because  cross-reference  information,  as 
with  the  wave  buoy  data  in  the  surface  wave  frequency  region,  does  not  exist.  These 
results  then  must  be  utilized  to  verify  the  ability  of  the  MRLS  algorithm  to  identify 
the  presence  of  low  frequency  energy  in  the  data  set.  In  this  respect  the  plots  show 
some  exciting  results.  Each  of  the  plots  shows  low  frequency  events  of  durations  up 
to  45  minutes.  These  events  track  in  frequency  and  the  dynamics  are  clear.  The  solid 
arrows  in  Fig  5.4  point  out  just  a  few  of  these  events. 

In  each  plot  there  exist  traces  that  show  large  jumps  in  spectral  level  and  have 
a  completely  different  character  than  the  other  spectra,  the  shaded  arrows  in  Fig  5.4 
highlight  these  traces.  The  character  is  induced  by  the  large  track  transitions  noted 
in  the  previous  section.  These  transitions  appear  as  steps  in  the  decimated  time 
series  and  thus  induce  ringing  in  the  spectra  when  encountered  in  the  processing.  The 
algorithm  quickly  adapts  to  the  new  level  and  continues  with  useful  spectral  estimates 
after  the  steps.  The  validity  of  spectral  estimates  in  the  vicinity  of  these  steps  must, 
however,  be  held  suspect.  None  the  less,  the  algorithm  allows  the  maximum  amount 
of  information  to  be  extracted  from  the  data  set  despite  the  contamination. 
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Frequency  (X  0.01  Hz) 

Figure  5.4:  Dynamic  spectral  waterfall  display  at  internal  wave  frequencies 
for  arrival  A  station  J. 


Figure  5.5;  Dynamic  spectral  waterfall  display  at  internal  wave  frequencies 
for  arrival  B  station  J. 
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Figure  5.6:  Dynamic  spectral  waterfall  display  at  internal  wave  frequen- 
fCiesfor  arrival  C  station  J. 
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Some  of  the  structure  of  the  arrival  A  and  C  dynamic  spectra  could  be  explained 
away  by  the  track  contamination  described  earlier.  However,  arrival  B  has  good  track 
quality  and  the  short  duration  events  are  evident  in  this  track.  This  would  imply  some 
physical  significance  but  the  source  is  not  necessarily  the  effects  of  internal  waves  (i.e. 
interference  of  a  closely  spaced  arrival  might  be  biasing  the  arrival  time  estimate  in 
some  long  term  measurable  manner).  The  ability  of  the  MRLS  method  to  identify 
low  frequency  dynamics  is  evident.  This  is  a  significant  result.  All  three  arrivals 
also  show  in  the  plots  energy  of  very  low  frequency  near  DC.  This  is  investigated  by 
applying  a  second  Chebychev  low  pass  filter  with  a  corner  frequency  of  0.0018  Hz  to 
the  arrival  series  already  decimated  by  a  factor  of  10.  These  series  are  decimated  by 
a  factor  of  10  again  yielding  a  time  series  for  very  low  frequencies  of  approximately 
100  points.  The  MRLS  algorithm  is  applied  to  this  series  with  a  forgetting  factor  of 
one  but  all  other  parameters  as  in  the  first  decimation  case. 

Figure  5.7  indicates  the  results  of  this  processing.  All  three  arrivals  show  peaks 
very  near  DC.  Arrivals  A  and  C  show  two  peaks  of  higher  frequency.  These  may 
associated  with  the  low  frequency  effects  of  the  track  transitions  as  discussed.  The 
lowest  frequency  peak  is  likely  attributed  to  physical  phenomena.  The  utility  of 
the  MRLS  technique  developed  has  been  demonstrated.  It  is  useful  to  note  that 
this  algorithm  can  also  use  complex  data  without  modification.  The  data  set  has 
phase  information  which  could  be  used  directly  in  the  spectral  estimation  if  a  reliable 
method  of  estimating  the  phase  at  the  LMS  filtered  arrival  time  could  be  determined. 
This  would  provide  a  significant  improvement  in  the  frequency  resolution  capability 
on  this  data  set. 
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C.  RECOMMENDATIONS 


The  following  items  are  suggestions  that  arise  as  result  of  working  with  the 
MBTE  data  and  would  assist  in  data  interpretation  if  implemented. 

1.  Calibrate  the  recorded  data  to  actual  sound  pressure  levels  received. 

2.  Reprocess  the  data  to  make  use  of  more  dynamic  range  on  the  A/D  conversion. 
The  present  integer  values  in  the  data  do  not  exceed  1000  and  are  more  often 
in  the  range  of  200  to  300.  This  corresponds  to  using  approximately  12%  of  the 
available  dynamic  range  of  a  12  bit  A/D  converter. 

3.  Perform  a  simple  ambient  noise  analysis  of  each  of  the  stations  to  determine 
the  ambient  noise  environment  that  the  acoustic  receivers  are  operating  in. 

4.  Track  the  dynamics  of  the  acoustic  carrier  frequency  in  the  raw  acoustic  data 
of  each  station  to  determine  carrier  stability.  This  is  easily  accomphshed  as 
part  of  the  ambient  noise  analysis.  It  would  have  been  useful  to  have  recorded 
the  actual  output  of  the  transmitter  from  a  receiver  placed  close  by,  that  is  to 
monitor  any  drifts  in  its  performance. 

5.  Compare  the  results  of  conventional  matched  filtering  with  the  Hadamard  trans¬ 
form  matched  filtering  to  see  if  better  resolution  on  the  arrival  peaks  can  be 
achieved.  The  better  the  resolution  on  the  arrival  peaks  the  better  any  tracker 
will  perform  on  the  data  set. 


With  a  view  to  future  research  based  on  the  results  of  this  thesis  work,  there 
are  some  difficulties  with  the  separate  LMS  and  MRLS  techniques  that  should  be 
investigated.  The  primary  difficulty  is  with  initialization  of  the  LMS  routine  and 
then  control  of  which  sub-arrival  that  the  algorithm  chooses.  Note,  both  the  LMS  and 
MRLS  techniques  utilize  a  set  of  tap  weight  coefficients  for  an  AR  process.  The  MRLS 
technique  is  a  high  resolution  principle  component  technique  and  the  LMS  is  a  fast 
efficient  adaptive  technique.  Both  of  these  may  be  combined  through  the  tap-weight 
vector  to  produce  a  single  implementation  which  could  prove  to  be  the  phase-lock 
loop  of  peak  tracking.  This  could  be  effected  by  selecting  a  peak  from  a  peak  set  that 
maximizes  a  principle  component  in  the  covariance  matrix  of  the  MRLS  technique 
and  then  adapts  the  tap-weights  of  the  predictor  to  track  the  component.  This 
approach  should  produce  an  algorithm  that  is  very  responsive,  efficient,  controllable 
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and  optimal  in  that  maximizes  some  selected  criterion.  Results  of  research  into  this 
aspect  of  the  signal  processing  could  produce  some  impressive  results. 
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VI.  CONCLUSIONS 


A  goal  of  this  thesis  work  was  to  develop  algorithms  to  improve  the  arrival  time 
tracking  of  the  Hadamard  matched  filter  output  of  the  MBTE  data.  In  addition, 
dynamic  spectral  estimation  of  the  possible  nonstationary  arrival  tracks,  in  the  surface 
and  internal  wave  frequency  regions  was  desirable.  These  efforts  have  yielded  the 
following  results: 

1.  Multipath  interference  was  shown  to  exists  in  the  station  J  arrivals  processed 
for  the  MBTE  data  set.  These  closely  spaced  arrivals,  in  most  instances,  are 
pajtially  resolved.  Their  presence  is  also  predicted  by  3-D  modeling  [Ref.  4]. 

2.  The  anomalous  arrival  fluctuation  levels  noted  in  [Ref.  1]  are  a  result  of  this 
multipath  interference. 

3.  Working  with  the  assumption  that  the  arrival  paths  are  fairly  stable  and  that 
close  phase  relationships  cause  high  amplitude  fluctuations  between  the  arrivals, 
an  LMS  tracking  technique,  independent  of  peak  amplitude,  was  developed  and 
implemented  showing  superior  performance  over  the  original  arrival  tracking 
technique  for  station  J  arrivals.  This  algorithm  minimizes  the  multipath  effects. 

4.  The  presence  of  the  surface  wave  component  was  verified  by  conventional  spec¬ 
tral  processing  of  the  LMS  arrival  track  error  signals  for  station  J. 

5.  A  high  resolution  MRLS  spectral  estimation  technique,  specifically  developed 
to  process  nonstationary  arrival  tracks  for  this  data  set,  revealed  low  frequency 
periodic  energy,  in  the  internal  wave  spectral  region,  in  all  three  arrivals  from 
station  J.  This  energy  cannot  be  attributed  to  internal  waves,  but  the  algo¬ 
rithm  demonstrates  the  capability  of  estimating  the  dynamic  spectra  of  these 
processes. 
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APPENDIX  A:  LMS  COMPUTER  CODE 


The  routines  in  this  appendix  implement  the  LMS  tracking  algorithm  in  the 
MATLAB  program  language  and  its  MEX  file  FORTRAN  interface.  Briefly  these 
routines  are: 

1.  PROGRAM  TLMSI  -  The  interactive  input  for  the  main  tracking  program 
TLMS. 

2.  PROGRAM  TLMS  -  The  main  line  LMS  tracking  code. 

3.  PROGRAM  SFRM  -  A  support  routine  for  TLMS.  The  2nd  derivative  compu¬ 
tation  for  locating  local  maxima  in  the  track  window. 

4.  PROGRAM  RESHAPE  -  A  support  routine  for  TLMS.  Rearranges  the  number 
of  rows  and  columns  of  a  matix  into  the  desired  dimensions. 

5.  PROGRAM  GETFRM  -  A  FORTRAN  routine  written  for  the  MATLAB  MEX 
file  FORTRAN  interface  designed  to  read  data  fron  the  Bernoulli  Box  mass 
storage  disk  drives  that  contain  the  data. 

6.  PROGRAM  GETFRMG  -  A  FORTRAN  routine  that  implements  the  MEX  file 
interface  between  the  FORTRAN  subroutine  GETFRM  and  MATLAB. 


A.  PROGRAM  TLMSI 

X  lot«r>cti*«  input  rootin*  for  th«  LHS  trnckor  TUfS.  Inpnto  nro  oolf 
X  axplanltory  and  aaplifiad  in  tha  TLIIS  rontina. 

X  Craatad  by  :  Ltd)  E.K.  Cbaulk  17  Qctobar  1990  (axactly  1  yaar  aftar) 

X  Tbaais  (tba  big  Qoaka) 

diary  a:\jl418d.aac;  diary  off; 

ib8*inpnt(’Entar  Start  Bin  nnabar  for  tha  tracking  aindoa  (1--134}  ■  '); 
ibu- input (’Entar  End  Bin  lumbar  for  tba  tracking  aindoa  Oibs-134)  ■■  ’); 
itpainpnt( ’Entar  FFT  Intarpolation  langtb  in  track  nindoa  (&13  ate)  ■  ’); 
ordainpntCSalact  LHS  filtar  ordar  {>64  prafarrad)  «  ’); 

■ntainpntCSalact  tba  tima  filtar  adaptation  parametar  mn  C'.OOl)  «  ’); 

■naa input (’Salact  tba  amplitnda  filtar  adaptation  paramatar  mo  C'lE-8)  >  '); 
tfxainpntC’Salact  tba  Lock  tima  for  initialization  >  ’); 
pltaO  X  Input  not  naad 

X  pltainputC ’Salact  tha  procasaad  data  output  moda  •  ’); 
ibufainputC 'Entar  tha  numbar  of  framaa  to  buffar 
filiainpot( 'Entar  tba  input  fila  nama  for  procaaaing  ■  '); 

X  Form  input  array 

xd«[iba  ibn  itp  ord  not  mua  tfx  pit  ibuf] ’ ; 

X  Start  tba  trackar 
tlmx(xd,f  ill)  ; 
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B.  PROGRAM  TLMS 

fanctioB  tluCzd.f ill) 

X  Thia  roatla*  ijq>l«Mnta  a  laaat  Baan  aqaaxaa  algorltba  for  arrival 
X  tracking  of  data  f ron  tha  Nontaraj  Ba;  Toaograpkj  Ezpariaaat .  Tha  input 
X  array  10  and  tha  fila  nana  FILI  ara  inpnta  fron  tha  input  routina  TLMSI 
X  ahich  ia  intaractiva.  Thia  iaplasantation  alloaa  a  aingla 
X  paraaatar  in  10  to  ha  aodifiad  uithout  having  to  go  through  tha  antira 
X  intaractiva  input  routina. 

X  Craatad  by  :  Ltd)  E.K.  Chaulk  17  Octobar  1990  (azactly  1  yaar  aftar) 

X  Thaaia  (tha  big  Quaka) 

X  Initializa  variablaa  for  tha  algor ithn 
iba>zd(l};  X  Starting  bin  for  track  vindoa 

ibn>zd(3);  X  bin  for  tha  track  uindoa  (  t  pointa  a  pouar  of  2!!) 
itp>xd(3};  X  FFT  intarpolation  langth  (poaar  of  3  io  S12) 

ord>zd(4) ;  X  Salactad  filtar  ordar 

■at>zd(S>;  X  idaptation  paranatar  for  tha  arrival  tiaa 

■ua>zd(6) ;  X  idaptation  paranatar  for  tha  amplituda 

tfzazd(7);  X  Sat  tha  avaraga  arrival  tint  of  intarast 

plt>xd(S);  X  lot  inplanantad  (to  ba  unad  for  dabngging  coda) 

ibufvzdO);  X  lunbar  of  franaa  to  buffar  from  disk  (nachina  nanory  dap.) 

inataO;  X  Valua  for  UIS  init  of  nav  arrival  tiaa  coaff  during  startup 

inaa>0;  X  Valua  for  LMS  init  of  nav  arrival  aaplituda  coaff  during  startup 

ichn>33;  X  Fortran  fila  I/O  channal  nuabar  for  fila  intarfaca 

ilan«134;  X  luabar  of  points  in  a  data  fraaa 

tpltaO;  X  Dagugging  paraaatar  for  routina  SFEM 

swO;  X  Initialization  conntar 

niniaord-1;  X  luabar  of  coaff s  raquiring  init 

nvibn'ibv^l ;  X  luabar  of  points  in  track  vindos 

inaB^az(aiza(aba(fili)})vi;  X  luabar  of  charactars  ia  fila  naaa 

Xlnitializa  output  variablaa.  Thasa  ara  tha  output  shift  ragistars 

tapvzaros (ord , 1 ) ; 

aat>vzaros(ord,l) ; 

p^>vzaros  (ord ,  1 ) ; 

att>zaros(ord,l) ; 

aaavzaros(ord,l} ; 

t^fvzaros(ord,l} ; 

aBpfazaros(ord,l) ; 

crssvzaros(ord,l) ; 

X  Mora  initializations 

otainat;  X  Init  first  arrival  tiaa  saight 

saainaa;  X  Init  first  aaplituda  saight 

yptaO;  X  Pradictad  arrival  init 

ypa«0;  X  Pradictad  a^lituda  init 

icataO;  X  Conntar 

itotaO;  X  Countar 

X  Main  tracking  loop 
shila  itot  <  13000, 

[zb  xp]agatfrB(lchn,f ili.inaa.ilan.ibuf) ;  X  Eaad  data  fraaas  froa  fila 

if  itot  av  0,  inaaa-inaB;  and;  X  Sat  val  for  raad  routina  aftar  first  raad 

X  Sat  track  vindos  and  arranga  buffarad  data 
ZBazB( ibs : ibn , ; ) ; 
zpazp( ibs : ibn , ; } ; 

ZBarashapa(zB,l ,a*ibuf) ; 
zpaTashapa(zp,l ,n*ibuf ) ; 

X  foliosing  cods  usad  during  initialization 
if  itot  <  ord. 
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for  i*l:lbQf, 
ltot>itot+l ; 

CtB,am,pa]Bofz»(zB,xp,B,i,itp,ibs  ,tplt} ;  X  Find  paaJui  2nd  dorr  rontin* 
if  itot  >  1. 

yptowt’*t^(itot-os:itot-l};  X  Prodict  with  arailoblo  saights  tiaa 
jpa>sa'*aq>(itat-ss; ltot-1) :  X  Saaa  for  aaplltuda 
and; 

if  itot  <  nini, 

[•t,itt]>Bin(ab8(tB-tfz));  X  Find  closaat  paak  during  inlt  tiaa 

alsa, 

[nt.ltt]>Bin(aba(tB-7pt)};  X  Find  cloaat  paak  aftar  init  tiaa 
and; 

[aa,lta]>Bin(aba^aB-ppa)) ;  X  l^litnda  closaat  paak 

X  Coapnta  arrors  for  aaight  npdata 
at»tM(itt)-ypt ; 
ta*aB(itt)-jpa; 

X  Updata  output  sbift  ragiatars 
t^(itat)>tB(itt) ; 
aaip(itot)>aB(itt) ; 
pap(itat}-pB(itt} ; 
att(ltot}aat; 
saa(ltot)>aa; 
tapf (ltat}*ppt ; 
a^if  (itot)*jpa; 
cras(itot)*icnt/ itot *100; 

X  Updata  aalghts  and  aztand  saighta  during  inlt 
if  ita  itt,  icntaicnt'^l ;  and; 
if  itot  >  1, 

ataat+(Bntaat) .*tBp(itot-au;itot-l) ; 
oaaaa-»(aniaaaa)  .aaap(itot-ss;itot-l) ; 
at* [at  ;  inat] ; 

aaafaa  ;  inaaj  ; 
and; 

saaBazCalzaCat)} ; 

[itt  tadtt)  ypt  at  lent/itot  Itotj  X  Output  displays  for  info 
[its  aaCitt)  ypa  aa  pa(itt)  i] 
and; 

X  Saaa  rasulta  in  diary  fila 

if  (ord-itot)  <  Ibuf,  lbof>ord-ltot ,  and; 
if  (ord-itot)  "■  0  , 

diary  on;  diapCCtapf  tq>  att  aapf  aap  aaa  pap  eras]);  diary  off; 
aa^iT(aiia(at)) ; 
ibnf«zd(9}; 
and; 

alsa 

X  Loop  uaad  aftar  initialization 
for  1*1: ibuf, 
ltot*ltot*l ; 

[ta,aa,pa]*sfra(za,zp,n,i,itp,ibs,tplt) ;  XGat  paak  locations  2nd  daria 

X  Pradictions  and  find  closaat  paaks 
ypt*st ’atapCord-ao^l :ord) ; 
ypa*Ba’*aap(ard-ar<'l  :ord)  ; 

[at ,itt]aain(abs(ta-ypt)) ; 

[aa,ita]>ain(abs(aa-ypa)) ; 

id»itt: 

if  itt  -•  ita, 
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X  Thasa  linaa  allow  aMplitnda  biaa  if  dasirad. 

[zdwJB  ldaM3amis(ab8([tB(itt)-tBpf(ord)  tB(ita)-tBpf  (ord)3)}  : 

X  if  idoa  2,  id'ita;  alsa,  id>itt;  and; 

alsa 

lcnt>icntal ; 
and; 

X  Coapnta  arror  for  naight  npdataa 
at"t»(id)-ypt ; 
aa»a»(id)-ypa; 
j*raa(i,ord) ; 

X  Updata  ontpnt  shift  ragistars 
tBp>[tap(2:ord)  ;  tn(itt)] ; 
aMp>[aap(2:ord)  ;  aa(ltt)] ; 

Fap*[pap(2:ord)  ;  paCitt)]; 
att»[att(2:ord)  ;  at]; 
aaa«[aaa(2:ord)  ;  aa] ; 
tapf»[tapf(2;ord)  ;  ypi^ : 
aBpfa[aapf(2:ord}  ;  n>*3 : 
craa>[crsB(2:ord)  ;  icnt/itot*100] ; 

X  1X3  tap-waight  npdataa  for  both  aaplitnda  and  tiaa 
BtastaCantaat)  .•tsq>(ord-sw+i  ;ord)  ; 
wa>wa+(ma*aa)  .aaapCord-aw^l  ;ord) ; 

[Itot  icnt/itot*100  tfz  jpt  tn(itt)  at]  X  Scraan  info  npdataa 

if  raa(itot,ord)  ■■  0, 

X  Sa*a  data 

diary  on;  diapCCtnpf  tap  att  aapf  aap  aaa  pap  eras]);  diary  off; 
and; 
and; 
and; 
and; 

X  Claan  np  aftar  last  data  raad 
iandaord-raaCltot ,ord)al ; 

diary  on;  diap([tapf (iand;ord}  t^>(iand:ord)  ... 
att ( iand : ord)  aapf ( iand : ord)  aapCiandiord)  aaa(iaad:ord)  ... 
pi^(iand;ord)  craa(iand:ord}]) ;  diary  off; 

and; 


C.  PROGRAM  SFRM 

fnnetion  [tB,aB,pa]>sfrB(zB,zp,ln,fn,innB,lt ,bt ,tplt) 

X  CTM,l]l,Plf]>SFUI(XH,ZP,LI,nr,lT,BT,n>l.T)  This  fnnetion  parforas 
X  diffarantiation  on  tha  intarpolatad  arrival  to  datoraina  tha  local  Haziaa  and 
X  ratnms  tha  list  of  candidata  arrival  tiaas  with  tha  corrasponding  aaplitnda. 
X  U  and  IP  ara  tha  aagnitnda  and  phasa  data  arrays.  LI  is  tha  fraaa  langth, 

X  PI  is  tha  fraaa  nnabar  in  tha  prasant  bnffar,  IIUH  if  tha  fraaa  nnabar  for 
X  display  in  dabngging,  IT  is  tha  FFT  Intarpolation  langth  and  BT  is  tha 
X  basa  tiaa  point  nnabar  for  convarsion  to  absolnta  tias  dalay.  TPLT  is 
X  tha  dabngging  and  plotting  inpnt .  Tha  ontpnts  TH,  iH  and  PH  ara  tha  paah 
X  arrival  tiaa,  arrival  aaplitnda  and  arrival  phasa  raspactlvaly . 

X  Craatad  by  ;  Ltd)  E.l.  Chanlk  17  Octobar  1990  (azactly  1  yaar  aftar) 

X  Thasis  (tha  big  Qnaka) 

X  Usa  data  fraaa  trach  window  for  aagnitnda  and  phasa  data 
vaza((fn-l)*ln+l :fn*ln) ’ ; 
vp»zp{  (fn-Dvln+l  :fn*ln)  ’ ; 
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X  Istatpolota  both  Bagnitndo  and  pbaa*  valnaa 
x»f  iiitr(»,lt) ; 
zp^intr(Tp,it) ; 

X  Uaa  2iid  darlTativa  to  locata  paaka  and  throa  aaa;  local  Bininaa 

*1>[0  ;  -Bln(0,diff(aign(diff (z))))  :  0  ]; 

tB>f ind(Tl>0) ; 

aB>z(ta) ; 

pa>zp(tB} ; 

X  Tbia  aaction  of  coda  could  ba  oaad  to  inplaaant  aaplitnda  biaaing  of  tha 

X  paaka 

Xnallaaaz ( a iza ( am) ) : 

Xittt>f ind(aB  >  0.7B*anB(aB)/nall) ; 

XaB'aadttt)  ; 

Xl«“pB(ittt)  ; 

XtB«tB(ittt); 

X  Dabugging  Coda  fox  plotting  routina  raaulta 
if  tplt  >  0, 

t-((bt+((0;lii-l)))  •(1.9375/124))’ ; 
tl-C(bt+(0;it-l)  ■•(Wit))  .•(1.9375/124))’; 

•lazaroa(Tl) ; 

•1 (tB)aonaa(tB) ; 
if  tplt  —  1, 

plot(t.T,’ ; ’ .tl.z. ’-’ ,tl ,(t1 .•Baz(»)). ’-») 

tart>aprintf ( ’Intarpolatad  Irrivala  sitb  paaka  Fraaa  «  Xg’.inua); 

alaa 

plot(tl ,z, ’-’) 

tazt>aprintf( ’Intarpolatad  Irriaala  Fraaa  ■  Xg’.inna); 
and; 

zlabaK’Tiaa  Oalaj  (aac)’) 
jlabaK  ’laplitnda’) 
titla(tazt) 
and; 

tBa(bt^(tir-l)  .aCln/it))  .*(1.9373/124) ;  X  Conaart  to  actual  tiaa  dalay 
and; 


D.  PROGRAM  RESHAPE 

function  y  >  roabapa(z ,a,n) 

X 

XY>R£SUPE(Z,N,I)  ratuma  tba  H-by-l  aatriz  aboaa  alaaanta 
X  ara  takan  coluanaiaa  from  1 .  In  arror  raaulta  i*  I  doaa 
X  not  haaa  Hal  alaaanta.  Tbia  la  a  built  in 
X  lUTLlB  function  oaad  to  raaiza  aatricaa 

LK,an]  >  Blza(z) ; 

if  BB*nn  B*n 

arror( ’Matrix  Boat  bava  H*l  alaaanta.’) 
and 

y  a  zaroata.n); 
y(:)  -  z; 

E.  PROGRAM  GETFRM  AND  GETFRMG 

c 

C  Tbia  ia  a  FOETRU  anbrontina  naad  to  opan  a  data  fila  and 
C  raad  a  nnabar  of  pointa.  Tbia  routina  ia  uaad  by  tha  MEI  fila  intarfaca 
C  uaa  of  lUTLlB  and  ia  baaically  af  atandard  FORTRII  routina.  In  ordar 
C  for  it  to  Bork,  It  auat  ba  coapilad  with  tha  KEl  conand  and 
C  HEX  librariaa  anppliad  nith  HITLIB  an  Ball  aa  an  intarfaca  routina, 

C  in  tbia  caaa  GERFRMQ.FOR.  Tba  ouputa  froa  froa  tbia  routina  ara  TH 
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C  shlcli  la  tha  magnitnd*  array  raad  in  Iroa  tha  Ilia,  TP  ahich  ia  tba 
C  phasa  array  raad  in  from  tha  fila  and  an  arror  paramatar  EZ.  Tha  inpnta 
C  ara  tha  fila  FO&TKII  channal  ntuabar  IP,  Tha  fila  na>a  FP,  tha  anmbar 

C  of  charactara  in  tha  fila  naaa,  Tha  langth  of  a  data  fraaa  L  and  tha 

C  noBbar  of  fraaaa  to  road  in  I.  Paraaatar  paaaing  throngh  tha  intarfaca 
C  rontina  la  tricky  so  plaaaa  raad  tha  KEZ  fila  Inforaation  in  HATLIB 
C  for  dataila.  Iota  if  tha  nimbar  of  charactars  in  tha  fila  naaa  is 

C  positiaa  tha  fila  is  opan,  is  it  is  nagatiaa  tha  fila  ia  accsssad 

C  and  if  It  is  0  tha  fila  ia  closad. 

C 

C  Craatad  by  ;  Ltd)  E.l.  Chanlk  17  Qctobar  1990  (azactly  1  yaar  aftar) 

C  Thasis  (tha  big  Qnaka) 

SUBROUTIIE  GETFEMCYK,  TP.  EZ.  IP,  FP,  DP.  L.  ■) 

»ElLa8  TJICI).  TP(1),  EZ.  IP,  FP(1).  DP.  L,  ■ 

IITEGEKa4  HOUS.NIIUTE.IOS .lUIIT 
IITEGER*2  UG, PHASE 

CHAUCTER*60  HEIDRC 
CBlRlCTERalS  H1C.H2C 
CBiRiCTER*eO  FHIiM 

lUIIT-HTdP) 

IF(IIT(OP).Eq.O)THEI 

CLOSE(IUIIT) 

RETUU 

EIDIF 


IF{IIT(DP).GT.O)THE* 

DO  10.I-1,IIT(DP) 

FILIi)((I;I)aCBiR(IfT(FP(I))) 

10  COITirUE 

OPEIdUIIT.FILEaFILIU.STlTUSa'aLD’  ,ERR-999) 

REiDdUIIT ,  >  (160)  >  ,ERR»999)BE1DRC 

RElDdUlIT  ,*(116,13,19,13)*  ,ERRa999)HlC  .HOUR  .B3C  .KIIUTE 

HRITE(a,a)BElDRC 

HRITE(a,a)BlC, HOUR, H3C, KIIUTE 

VRITE(a,a) 

EIDIF 

HRITE(a,a}LaI 

DO  110  I>l,IITa*I) 

REiDdUIIT,  •)H1G.PB1SE 

HRITE(*,*}IUG,PB1SE 

TK(I)>IUa 

TPd)-PHlSE*3 . 141693E-3 
110  (niTIIUE 
EZ>0 
RETURI 

999  BRITE(a,a)*FUa  I/O  Error* 

EZ-1 

RETURI 

EID 
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PKOGRIH  GETFBIIG 

C  This  rontiit*  Bodifiad  froa  th«  axaq)!*  proaidad  in  tha  M£X  fila 
C  ijttarfaca  to  aoxk  aith  tha  GETFRM.FOR  anbrontlna.  Saa  tha  KATLiB 
C  dociiBaiitat ion  for  dataila  on  this  fila. 

C 

C  Editad  bj  :  Ltd)  E.I.  Chaolh 
C  Thasis 


tllCLUGE; ’FHEI.B’ 

IITERFACE  TO  SUBROUTIIE  GETFHM 

(YmCaslna],  YPPCaalna],  EZZCaalna], 

♦  IPCaalna]  ,  FPCaalna]  ,DP[»alua]  ,  LLCaalna],  ■■[anlue]) 
IITEGER*4  ym.  YPP,  EZZ,  IP.  FP,  DP,  LL,  II 
EID 

SUBROUTIIE  USRFCI 
+  [c .alias : ’usrfcn’] 

a  (ILHS.  PLBS[rafaranca] ,  IRHS,  P?RS[rafaranca]) 

IITEGERaZ  ILHS,  IRHS 
IITEGER*4  PLBS(*),  PRHSC*) 

IITEGERa4  CRTHAT,  REALP,  IHAGP,  GETGLO,  ALREAL,  ALIIT 
IITEGERaZ  BPCBK 
REALae  GETSCA 

IHTEGER*4  THFI.  YPP,  DP,  FP,  IP,  LL,  II 
IITEGERaZ  H,  I,  HDD 
REALaS  XL,  Zl 

IF  (IRHS  .IE.  S)  THEI 

CALL  MEXERK 'GETFRM  raqniras  fiaa  input  argoBants*) 
ELSEIF  (ILHS  .IE.  3)  THEI 
CALL  REXERRC 'GETFRH  raqniras  too  ontpot  argnaant’) 

EIDIF 

XL  a  CETSCA(PRHS(4)) 

XI  a  GETSCA(PRBS(S)) 

H  a  IIT(IL) 

I  a  IIT(XI) 

HDDal 

PLHS(l)  a  CBTIUT(K,l,0) 

PLBS(2)  a  CRTIUT(H,l,0) 

PLBS(3)  a  CRTHATCHDD.KDD.O) 

Ym  a  REALP(PLBS(1)) 

YPP  a  REALP(PLBS(2)) 

EZZ  a  REALP(PLHS(3)) 

IP  -  R£ALP(PRHS(1)) 

FP  a  REALP(PRHS(2)) 

DP  a  REALP(PRHS(3)) 

LL  a  REALP(PRBS(4)) 

II  >  REALP(PRHS(5)) 

CAtJ.  GEIFRH(Yra,YPP,EZZ,IP,FP,DP,LL,ll) 

lETUU 

EID 
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APPENDIX  B:  MRLS  COMPUTER  CODE 


These  routines  are  used  to  implement  the  MRLS  spectral  estimation  technique 
in  the  MATLAB  programming  language.  Briefly  these  routines  are: 

1.  PROGRAM  MRLS  -  This  is  the  main  computational  routine  for  the  spectral 
estimation  technique. 

2.  PROGRAM  PRED  -  This  is  a  MRLS  support  routine  for  extending  a  time 
series  using  a  set  of  prediction  error  coefficients. 

3.  PROGRAM  PERIOD  -  A  routine  for  computing  a  zero  padded  periodgram  of 
a  time  series. 

4.  PROGRAM  ARPER  -  A  routine  for  computing  the  power  spectrum  of  a  set  of 
prediction  error  coefficients. 

5.  PROGRAM  MODCM  -  A  routine  for  generating  the  forward- backward  data 
arrangment  for  the  modified  covariance  matrix. 


A.  PROGRAM  MRLS 

fosctioa  Cf .xuct) 

X  [F,T,A>RLS(X,0BD£R,IV,FF,I1.IZ,PLT.IXT}  This  foAction  rstoms  th«  AR 
X  co«fflci«ftts  for  tho  tho  soloctod  aodol  ORDER  in  roctor  A  using  ths 
X  BATKII  aodifiod  corurisncs  forvsrd*hsck«srd  asthod.  A  Inclndss  s  constant 
X  cosfflclsnt  1  in  ths  first  position.  Vsctor  Y  rstnms  ths  spsctrsl 
X  sstiasts  coapntsd  using  ths  A  ssctor  for  ths  normslixsd  frsqnsncj 
X  points  in  vsctor  F.  X  is  ths  dots  array  and  ORDER  is  ths  sslsctsd  ordsr  of 
X  ths  corrslatlon  aatriz.  IV  is  a  row  vsctor  sntry  with  ths  nnabsr  of  ths 
X  sigsBvalnss  to  bs  nssd  In  ths  cosfficisnt  calculation.  FF  is  ths  forgstting 
X  factor  which  vast  bs  bstvssn  0  and  1 .  II  is  a  factor  that  vill  alios  this 
X  algorltha  to  work  szactlj  liks  ths  HFBLP  tschniqus.  Unlsss  its  uss  is 
X  undsrstood  froa  inspscting  ths  rout  ins  alvays,  sst  this  valus  to  0. 

X  IZ  is  ths  nuabsr  of  spsctral  points  to  bs  coaputsd  froa  ths  A  vsctor. 

X  This  valus  is  ussd  to  zsro  pad  ths  output  apsctrua. 

X  PIT  dsfinss  at  shat  intsrval  plots  of  ths  procsss  ars  dsslrsd.  For  szsapls 
X  sotting  this  valus  to  1  sill  plot  a  spsctrua  aftsr  sach  updats.  Sotting  it  to 
X  nsgativs  valus  vill  output  ths  tap'vsight  spsctra  at  tho  dssirsd  intsrval. 

X  IXT  toll  ths  prsdiction  routins  hoo  aanj  point  to  sztsnd  ths  ssriss  forward 
X  and  backward . 


X  Crsatsd  by  ;  Lt(i)  E.K.  Chaulk  17  Octobsr  1990  (szactly  1  ysar  aftsr) 

X  Thssis  Prograa  (ths  big  Quaks) 

k’aiazCsizsCz)) ;  X  luabsr  of  points  in  ths  data  array 

nw^azCsizsCnv)) ;  X  lusdisr  sigsnvaluss  to  uss  in  th  procsssing 

nn^unsi ; 

X  Ths  asan  can  bs  rscuxsivsly  rsaovsd  by  thsss  statsasnts  and  othsrs  in  ths 
X  Bain  loop 

XzBSan^san(z(l  :ordsr^an)) ; 
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lz(l  :ord*rHiii)«x(l  :ord«z-*'nn)  ./zM>n; 


7S>Bo<lai(x(l  lordar+nn)  .order)  ;  t  Create  the  aodilied  coeerience  date  natrix 
be[z(order^l  :ordar+iui)  z(l:iui)]‘;  X  Create  the  deaired  reeponae  eector 
theye’eb;  X  Coapate  theta 

yaaya'eya:  X  Coapnta  the  initial  correlation  utrlz 

•*zaros(th) ;  X  Zero  the  eeight  eector 

X  Main  update  loop 
lor  i»order+na:k-l , 

X  RecuraiTe  nean  renoval  atatenents 
X  x»oan»(i*xnean+x(i+l)) ./(i+1) ; 

X  x(i+l)«x(i+l)/x»ean; 

X  Update  the  correlation  and  theta  natricea  nain^  the  reenraiee  relationa 
yd“[lliplr{x(i-order*l ; i))  ;  x(i-order+2:i+l)J ; 
ya“lf . eya+yd  >  eyd ; 

th»ll.eth+yd’e[i(i'*-l)  ;  xCi-order+l )]  ; 

X  Spectral  conpntation 

il  pit  >  0  t  ren(i-order-nn+l .pit)  ««  0. 

[n  d  T]>a*d(ya);  X  Singular  Value  deconpoaition 

X  Create  the  tap-seight  eectro  Iron  aigenTalnaa  and  eigeneectora 
lor  j»l:n»».  »1( : . j)«T( : .neC j)) . /d(n»( j) .neC j)) ;  end; 
lor  j«l:nTT.  b«u+»1 ( : , j)» (t( : .n»( j) ) ’eth) ;  end; 
a»  [1 ;  -«]  ; 

il  nxt  >  0. 

z«pred(x(i-order-nn-plt+2 : i+1) .a.nxt) ;  X  Extend  the  data  eia  linear  pred 
■B>aaz(aize(z)) ; 

B^hanaingCnB)  > ;  X  tpplj  a  Baming  aindoB 

[1 .y] Bper iodCz . *b  .nz) ; 

diary  on;  diap(y’);  diary  oil;  X  Save  the  ontput 

X  Plot  the  reaulta 
plotd.y) 

tezteaprintK ’Point  t  Xg,  nai  Xg’ .i+i  .Bax(ya)) ; 
title(text) 

ylabelC’ Magnitude  (dB)’) 

zlabeK ’Percentage  ol  Sanpling  Fra<ioency  (Hz)’) 
vezeroa(th) ; 

elae. 

X  Coapnta  the  tap-Baight  apectrun  inatead  ol  the  MRLS  apectrua 
Cl.y]aarper(a.nz) ; 
plotd.y. ’-.’) 

tezteaprintK ’Point  t  Xg’.!-*!); 
title(tazt) 

ylabelC ’Magnitude  CdB)’) 

ZlabeK ’Percentage  ol  Saapling  Frequency  (Hz)’) 
eezeroa(th) ; 
end; 

X  Sava  reaulta  il  deaired 
X  diary  on 

X  diap(a) 

X  diary  oil 

X  aeta  plotaa.pl 
and ; 
end ; 


B.  PROGRAM  PRED 

lunction  yeprad(z.a,n) 
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X  T>PEED(1,X,I)  Thia  fimctloix  is  ossd  to  aztsnd  ths  data  langth  of  a  saqaanca 
X  X,  basad  on  tba  coafflclanta  containad  in  tha  aactor  i.  I  aasplas 
I  ara  prapandad  to  tha  data  and  I  saaplas  ara  appandad  to  tha  data.  This 
t  rontina  is  dasignad  to  sork  sith  tha  pradiction  arror  coafficiants  prodncad 
X  txom  a  aariaty  of  prograns.  Tha  first  coafficiant  nnst  ba  a  coni.tant  valna. 

X  Craatad  by  :  Ltd)  E.X.  Chaolk  20  loTaabar  1983 

X  Thasis 

■■■az(siza(z)) ; 
l^az(siza(a))  : 
a— a(2:l)  ; 
af>flipad(a) ; 

1-1-1; 

y-[zaros(l ,n)  z  zaros(l.n)]; 
for  i»0:n-l, 

X  [i  n-i  n-i+1  n-i+1  i+«+l  n+i-l+l  n+i]  X  Indaz  check  for  debug 
y (n-i)»y(n-i+l :n-i+l)»a; 
y (i+n+li-y Cn+i-l+l :«+i)*af ; 
end ; 

end; 


C.  PROGRAM  PERIOD 

fonction  [f.y]  “  pariod(z,z) 

X  [F,Y]-PEKIOD(X,Z)  This  function  coapntea  tha  pariodograa  of  fft  sidth 
X  Z  of  tha  data  saqnanca  X.  Tha  noraalizad  fraqnancias  ara  ratnmad  in 
X  F  sith  tha  posar  spactral  astimata  in  aactor  y.in  dB's.  Hindoss  nnst  ba 
X  appliad  saparataly. 

X  Craatad  by  ;  Ltd)  E.X.  Chaolk  B  loaaabar  1989 

X  Thasis 

sb-2.0; 

if  z  <  0,  sa-1.0;  and; 
z-abs(z) ; 

n-2*(l*fiz(-.0000001+(log(»az(siza(z)))/log(2)))); 

I-0:(z/2-l); 
f-d./<z/2)./2]; 
z(z)-0.; 
y»fft<z) ; 
y-y(l:z/2); 
y-(s«/n) .eabsCy) ; 

Xy-y ./nazCy) ; 
y>uz(y,  .00000001)  ; 
y-20 .0 .eloglO(y) ; 

•nd ; 


D.  PROGRAM  ARPER 

function  [f.y]  -  arparCa.z) 

X  [F,T]-1PER(1,Z)  This  function  co^putas  tha  Posar  Spactral  Density  in  dB’s 
X  froa  the  sector  1  of  pradiction  arror  coafficiants  suppliad. 

X  Tha  salua  of  Z  indicatas  tha  fraquancy  rasolution  and  is  nsad 
X  to  zaro  pad  tha  output .  Tha  noraalizad  fraquancy  saloas  ara  ratnmad 
X  in  F. 

X  Craatad  by  ;  Ltd)  E.X.  Chanlk  20  losaabar  1989 

X  Thasis 

X 
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n«2-(l+f  ix{-. 0000001 •••Clog(Bai(8iz«(*)))/log(a)))); 
f-C(0;(2/2-l))./(2/a)./2J'; 

•««(l);a(l)>l; 

•(2)a0. ; 

j»B.*(2.0.«»b»(fft(*)) ./n) .-(-1) ; 

■J’mMxlj,  .00000001)  ; 

Xy*y ./(■*2{j)) ; 
y«10.0.»logl0(y) ; 

y«yCl:2/2): 

•ad; 


E.  PROGRAM  MODCM 

function  y^BodcaCz.ip) 

X  T>HOOCI((X<IP)  Tlila  program  simply  builds  a  matrix  in  th«  modified  coTariiuace 
X  data  arrangement  trom  a  rso  data  vector  input.  I  is  tbe  data  array  and  IP 
X  ia  the  order . 

X  Created  by  :  Lt(l)  E.X.  Chaulk  20  lovember  1989 

X  Thesis 

X 

n>max(8ize(x)) ; 
np«B-ip; 

yezeros(2enp,ip) ; 
for  i>l:np, 
kel : ip : 

y (i,k>»x(i+ip-k) ; 
y(i*np,k)“x(i*k) ; 
end 

end 
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APPENDIX  C:  LMS  TRACKER  OUTPUT 


This  appendix  contains  the  detailed  output  of  the  LMS  tracking  algorithm  for 
arrivals  A,  B  and  C  of  station  J  for  the  14  Dec  1988.  The  figures  are  arranged  in 
consecutive  groups  of  three  for  comparison  of  each  of  the  tracked  values  for  the  three 
arrivals.  The  seven  outputs  are; 

1.  The  LMS  filter  arrival  time  predicted  track. 

2.  The  measured  location  of  the  closest  arrival  time  peak  to  the  predicted  value. 

3.  The  difference  between  the  predicted  and  measured  arrival  times,  termed  the 
arrival  time  error  track. 

4.  The  LMS  predicted  amplitude  at  the  arrival  time  peaks. 

5.  The  measured  amplitude  at  the  arrival  time  peak  closest  to  the  LMS  predicted 
arrival  time  track. 

6.  The  measured  phase  at  the  arrival  time  peak  closest  to  the  LMS  predicted 
arrival  time  track. 
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Figure  C.l:  Arrival  A  LMS  predicted  time  track. 
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Figure  C.2:  Arrival  B  LMS  predicted  time  track. 
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Figure  C.3:  Arrival  C  LMS  predicted  time  track. 
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Figure  C.9:  Arrival  C  arrival  time  error  between  prediction  and  closest 
peak. 
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Figure  C.IO:  Arr'val  A  LMS  predicted  amplitude  track. 
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"ure  CM4:  Arrival  B  peak  amplitude  value  at  closest  peak  to  LMS 
dieted  arrival  time. 
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Figure  C.15:  Arrival  C  peak  amplitude  value  at 
predicted  arrival  time. 
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Figure  C.19:  Arrival  A  Phase  values  at  closest  peak  to  LMS  predicted 
arrival  time. 
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Figure  C'.20:  Arrival  B  Phase  values  at  closest  peak  to  LMS  predicted 
arrival  time. 
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Figure  C.21:  Arrival  C  Phase  values  at  closest  peak  to  LMS  predicted 
arriv'al  time. 
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APPENDIX  D:  MATCHED  FILTER 
CORRELOGRAMS 


These  are  the  16  correlograms  that  constitute  the  matched  filter  output  for 
Station  J,  14  Dec  88.  Included  are  LMS  and  low  pass  filtered  track  overlays  for 
the  A,  B  and  C  arrivals.  These  plots  validate  the  arrival  tracks.  Note,  the  original 
resolution  of  the  matched  filter  output  would  have  only  produced  124  pixels  for  the 
greyscale.  To  get  the  full  resolution  of  the  page  each  sequence  period,  or  trace,  was 
FFT  interpolated  from  124  points  to  512  points.  This  stretches  the  plots  in  the 
X-direction  and  allows  a  clear  picture  of  the  underlying  dynamics. 

Note,  the  last  9  to  10  minutes  of  Fig  D.16  show  a  non-zero  output  for  each  track 
since  the  algorithm  works  with  the  peaks  of  data  no  matter  how  small  they  are. 
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Figure  D.l:  Correlograin  7^1 
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Figure  D.3:  Correlograin  #3 
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Figure  D.6:  Correiogram 
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Figure  D.7;  Correlogram  #7 
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Figure  D.8:  Correlofrranj  #8 
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Figure  D.IO:  Correlogram  ^10 
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Figure  D.ll:  Corielogram  #11 
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Figure  D.12:  Correlogram  #12 
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